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FOREWORD 


The research program described herein was 
sponsored by the National Aeronautics and 
Space Administration, Lewis Research Center, 
Cleveland, Ohio, 44135, under contract 
NAS3- 11226. The study was conducted during 
the 14-month period beginning 22 May 1968 
and ending 22 July 1969. NASA Project 
Manager was Dr. R. J. Priem. The Rocketdyne 
program manager was Mr. T. A. Coultas, with 
Dr. C. L. Oberg as responsible engineer. 

Dr. Oberg was largely responsible for form- 
ulation of the variational-iterational 
technique; Dr. T. L. Wong was responsible 
for numerical solution of the resultant 
equations, and Dr. R. A. Schmeltzer devel- 
oped the secular equation method. 

This report has been assigned the Rocketdyne 
report No. R-8076. 
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ABSTRACT 


As a means of studying the effects of baffles on 
acoustic modes of combustion instability, analy- 
tical methods have been developed for calculating 
the wave motion in closed, baffled chambers with 
rigid or nonrigid boundaries. These methods 
encompass solving the wave equation for the 
baffled chamber by converting the differential 
equation and boundary conditions to an integral 
equation which, in turn, is solved by approximate 
means. A variational technique in combination 
with an iterated approximation was used to solve 
the integral equation. Numerical results were 
obtained for two-dimensional chambers containing 
one or several equal length and equally spaced 
baffles. The results show an essentially con- 
tinuous pressure distribution along the baffle 
tips. Requirements for continuity of velocity 
and energy flux are automatically met with this 
method. Furthermore, the effects of a single 
baffle on the stability of a chamber with non- 
rigid walls, i.e., gain/loss-type boundary condi- 
tions, has been successfully analyzed for one 
particular two-dimensional geometry. Thus, the 
ability to generalize the method for nonzero boundary 
conditions has also been demonstrated. Finally, 
convergence of the iteration scheme has not been 
proved but very good numerical results were ob- 
tained and it is likely that adequate convergence 
can be achieved. 
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NOMENCLATURE 


G a (r IV 


= coefficient defined by Eq. 12 (lb^/in?) 

= variational parameter 


= defined by Eq. 9 (lb^/in.) 

= sound velocity, in. /sec 


= Green's function for main chamber 


G bi/ r l r 0^ = Green ' s function for p th baffle compartment 


= o)/c, inch 

= (k 2 - q 2 ir 2 /w 2 ) 1/2 , inch ' 1 

= (k 2 - mV/W 2 ) 1/2 , inch ' 1 


= baffle length, inches 


= main chamber length, inches 


= (& + L) , total chamber length, inches 


= index, positive integral values 


= unit normal vector directed outward 
= pressure, lb^/in? 

= index, positive integral values 


position vector, inches 


= surface area, sq in. 
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CO 


total chamber- compartment interface area, sq in. 

t h. 

chamber- compartment interface area of p compartment, sq in. 

longitudinal velocity component, in. /sec 

3 

chamber volume, inch 
width of baffle compartment, inches 
width of main chamber, inches 
longitudinal position coordinate, inches 
transverse position coordinate, inches 

specific acoustic admittance of injector end (dimensionless) 
specific acoustic admittance of nozzle end (dimensionless) 

1 if v = 0, 2 if v / 0 

eigenvalue for reference volume, inch ^ 

{ J'i j>2 dV}/V, dimensionless 

baffle compartment index, positive integral values 
normal pressure gradient, lb^/in. 

3 

time averaged gas density, lb m /in. 
wW/c 

eigenfunctions for reference volume 
angular frequency, radians/sec 
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SUBSCRIPTS 


refers to main chamber 


b 

V 

0 


= refers to baffle compartments 
til 

= refers to y baffle compartment 


refers to source coordinates for Green's functions 


SUPERSCRIPTS 


s = 

overbar = 
(n) 

A = 


refers to source surface (interface) 

denotes particular index 
th 

refers to n iteration 

denotes maximum value retained in summation 
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SUMMARY 


Analytical methods have been developed for calculating the wave motion 
in closed, baffled chambers with rigid and nonrigid boundaries. Appli- 
cation of these methods to the design of injector-face baffles in liquid- 
propellant engines will provide significant insight into the effects of 
baffles on combustion stability. 

During this program, approximate solutions to the wave equation, with 
essentially continuous pressure distributions, were obtained for closed, 
two-dimensional chambers containing an unrestricted number of equal length 
and equally spaced baffles. Solutions were obtained by converting the 
wave equation and boundary conditions to an integral equation; this inte- 
gral equation was solved with a combination of variation and iteration 
methods. The mathematical techniques used to obtain these solutions apply 
equally well to cylindrical or annular chambers, and to unequal baffle 
lengths or spacing, although with some increase in complexity. As yet, 
however, numerical results have not been obtained for these more compli- 
cated configurations. The calculated frequencies and oscillatory pressure 
distributions for two-dimensional chambers are in good agreement with data 
obtained from bench-scale acoustic models. 

Calculations were also made for similar chambers with a loss- type boundary 
condition (admittance) at the nozzle end of the chamber, and a gain-type 
boundary condition at the injector end. These boundary conditions simulate 
pressure- coupled gains and losses. Results for the particular case con- 
sidered suggest that baffles improve stability in an engine by reducing 
velocity rather than pressure coupling. This indication is drawn from the 
fact that, in the pressure- coupled (acoustic- admittance) case considered, 
the addition of one baffle was found to worsen rather than improve stability. 
Nevertheless, the ability to employ the analytical approach with nonzero 
boundary conditions has been demonstrated. 

Convergence of the iteration scheme has not been proved, but very good 
numerical results were obtained. It appears likely that adequate convergence 
can be achieved. 
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INTRODUCTION 


Baffles are typically comprised of an array of blades, often straight, 
mounted on the injector face and projecting several inches downstream. 
Typical configurations may be found in Ref. 1 and 2. 

Currently, baffles are designed principally from consideration of the 
oscillatory velocity characteristics of observed or anticipated acoustic 
modes for the chamber. The baffle blades are generally positioned at or 
near the calculated location of the velocity antinodes (maxima); the 
antinodal locations are calculated from acoustic expressions that pertain 
to closed, unbaffled chambers. Thus located, the baffles are thought to 
"interfere" with the normal acoustic modes and thereby promote stability. 
Although this approach is often successfully applied, it is unsatisfactory 
because no consideration is given to the wave motion possible once the 
baffles are installed. 

Further, empirical design rules have been developed from motor firing and 
some acoustic modeling results, but these, too, are unsatisfactory because 
they lack quantitative character and are not always valid. Moreover, these 
results provide little insight into the manner in which stability is im- 
proved. Several of these techniques are described in Ref. 1 and 2. 

It is evident that a better fundamental understanding of the oscillatory 
characteristics of baffled chambers would contribute significantly to 
more nearly sound design techniques. Few previous analyses of these os- 
cillatory characteristics have been reported, perhaps because of the mathe- 
matical difficulties involved. One approximation used is obtained by first 
assuming the wave motion within the baffle compartments is only longitudi- 
nal, and then, matching the ratio of pressure and normal component of velo- 
city (admittance) at the interface between the compartments and the main 
chamber. This approximation yields a moderately good estimate of the 
frequency depression produced by baffles at the normal chamber frequencies. 
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However, in this approximation, none of the continuity requirements on 
pressure, velocity, and energy flux across the interface are met. Further, 
baffle spacing does not enter into the calculation. Therefore, this ap- 
proximation is lacking in many respects. 

Other attempts at analysis of the problem have been largely unsuccessful. 
Therefore, these provide little foundation on which to base another study. 
Fortunately, however, thorough analyses have been developed in the acous- 
tics literature for related kinds of problems. These do provide a satis- 
factory starting point for an analysis of baffled chambers. 

ANALYTICAL APPROACH 


The analytical approach is based on the assumption that the linear acoustic 
behavior of a baffled chamber adequately approximates the wave motion in a 
baffled combustion chamber. Although this approximation may be lacking in 
some respects, it has been effectively used for stability problems in the 
past, and further it is certainly an appropriate starting point for more 
thorough analysis. 

Thus, the effects of baffles were investigated by solving the wave equation 
for closed, baffled cavities. The effects of combustion driving and the 
various loss processes were simulated to some extent by employing gain/ 
loss-type boundary conditions. The effects of steady flow have been largely 
neglected; however, a modified ("virtual") boundary condition was employed, 
which partially accounts for this flow. 

Although the problem has been simplified to the extent that only the wave 
equation need be solved, the solution is still not straightforward. The 
presence of baffles (the boundary shape) precludes the use of the separation- 
of-variables technique to solve this equation. For similar problems Morse 
(Ref. 3, pg. 1039) suggests using an integral formulation of the problem 
and then solving the resultant integral equation by approximate means. 


4 



Accordingly, the wave equation and boundary conditions were rewritten as 
an integral equation* which was solved by two approximate means: an 

iteration-variation technique and a secular determinant technique. Sat- 
isfactory results were only obtained with the former method. 

The approximations were guided by, and tested against, frequency and 
oscillatory-pressure data from bench-scale acoustic modeling tests. It 
is clear from these data that the normal transverse modes of the chamber 
are distorted, but not eliminated, by the addition of baffles. Further, 
the baffle compartments and main chamber generally act as a closely coupled 
oscillatory system with a continuous pressure distribution near the inter- 
face between them. However, compartment modes, which are largely confined 
to the baffle compartments, also appear. These observations from modeling 
tests are compatible with engine test (hot-firing) data as well. 
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ANALYSIS 


INTEGRAL FORMULATION 


Because the wave equation cannot be solved by separation of variables for 
the baffled chamber, this equation and the accompanying boundary conditions 
were converted to an integral equation that was solved by approximate means. 

The steps involved in this conversion are described by Morse (Ref. 4, pagfe321); 
the Helmholtz equation, which is the wave equation for a harmonic time dependence 

V 2 P + k 2 p = 0 (1) 


may be rewritten as 


P(r) 






( 2 ) 


where G(r|rg) is a Green's function, which satisfies either the same boundary 
conditions as the pressure (p) , or a zero-gradient boundary condition. In 
addition. Green's function satisfies the differential equation 

V 2 G + k 2 G = -6(r - r Q ) ( 3 ) 

where 6(r - r^) is a Dirac delta function. 


The integral expression for pressure is used with separate Green's functions 
written for each baffle compartment and also for the main chamber. Each 
of these Green's functions is zero outside of the compartment to which it 
applies. However, the oscillatory pressure and normal component of velocity 
must be continuous across the interface between each region. Therefore, 
at this interface. 



*/ G a (? sl ? 0 J S(? 0 )dS 0 ■ -f V ? slV 5(? 0 )dS 0 = Pbp^s 5 

s s 

(4) 
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where G (r r^) is the Green's function for the main chamber, G. (r rj is 
a 1 0 ^ by 1 0 

the Green's function for the y baffle compartment, and 

- N-y p a (5 

(5) 

= -N-Vp b (^) 


Assuming the baffle compartments and main chamber have regular geometries, 
the necessary Green's functions can be readily developed. If the set of 
equations indicated by Eq. 4 are solved for £, then the pressure at any 
point within the combustion chamber can be calculated from Eq. 2. In addi- 
tion, a set of frequencies (eigenvalues) will be obtained from Eq. 4. 

Morse (Ref. 4, pg. 680) suggests approximate solution of Eq. 4 by a varia- 
tional technique; he also suggests a secular-equation technique for similar 
problems (Ref. 3, pg. 1040). From each of these methods, a characteristic 
equation is obtained which, in turn, must be solved for the eigenvalues or 
frequencies of the chamber. 

The variational-iterational technique will be described in some detail. A 
description of the secular-equation method will be given in Appendix A be- 
cause no satisfactory results were obtained with this method. 

VARIATION- ITERATION TECHNIQUE 

Morse (Ref. 4, pg. 680) developed a variational function from Eq. 4, al- 
though he was considering only two coupled cavities. An approximate 
expression for £, or rather the normal velocity, containing an unspecified 
parameter was inserted into the variational function and the variational 
procedure was used to optimize the parameter, i.e., select the "best" value. 
From this, a second-order estimate of the eigenvalue (frequency) and a 
first-order estimate of the eigenfunction (pressure) was obtained. 
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If this procedure is applied to the baffle problem with an approximate 
function of the form E, = A £ being used, where A is the variational 
parameter but E, is a still unspecified function of position, the vari- 
ational procedure leads to the result 

^(P a - P b ) ds = 0 (6) 

where 

Pa = ~/ G a ( ^s^0 ) T(? 0 } ^ 

s 

p b = v as 

s 

The mathematical . details concerning the development of Eq. 6 are summarized 
in Appendix A. Because E, is proportional to the normal component of velo- 
city, this procedure shows that the best approximate eigenvalue (frequency) 
is obtained by requiring continuity of the energy flux at the interface 
between the main chamber and the baffle compartments. Continuity of the 
normal velocity is automatically satisfied but pressure is not, unless the 
exact solution is used. 

During the initial part of this program, several approximations were tried 
for E,. (The overbar on E, will be subsequently omitted because the distinc- 
tion is unimportant.) However, in each case, the pressure match across 
the interface was poor. Consequently, an iteration procedure was devised 
to improve these initial estimates. 

The iteration was done with the aid of orthogonality properties. Green's 
functions were obtained as expansions in terms of orthogonal functions, 
and similar kinds of expansions were obtained for the pressures. The 
iteration was done by: (1) initially assuming a pressure distribution 
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on the baffle- compartment side of the interface; (2) calculating the cor- 
responding normal gradient (£) from the expansion ;(3) then calculating 
the pressure on the chamber side from that expansion for £; and (4) finally 
equating expansions for p a and p^ to obtain a new set of coefficients for 
p^ and complete the cycle. 

Convergence of this iteration scheme was not proved; however, good numerical 
results were obtained. 

GAIN/ LOSS- TYPE BOUNDARY CONDITIONS 

Nonzero, admittance-type boundary conditions can be added at each end of 
the chamber without difficulty. By defining Green's functions, which satisfy 
the same boundary conditions at the closed ends of the chamber as those 
satisfied by the pressure, the .foregoing equations may be used without 
change. 
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ANALYTICAL RESULTS 


TWO-DIMENSIONAL CHAMBER WITH RIGID WALLS 

With the variational method, continuity of the normal component of velocity 
and energy flux across the interface between the main, chamber and baffle 
compartments is automatically satisfied. However, continuity of pressure 
is not satisfied unless the exact solution is found. This interface, as 
sketched below, is, of course, simply a convenient surface on which to 
match solutions for the compartment and main chamber. 



In an effort to satisfy the pressure condition as well as possible, three 
different kinds of approximate functions were used. In the first of these, 
the pressure distribution on the main- chamber side of the interface was 
assumed to be of the form cos mny/W. Secondly, a similar distribution was 
assumed for the compartment side. Neither of these resulted in a good 
pressure match; therefore, the iteration scheme was developed to improve 
the agreement. 

MAIN-CHAMBER APPROXIMATION 

With the main-chamber approximation (i.e., assumed pressure profile on 
chamber side) the following characteristic equation was obtained (see 
Appendix B) for equal length and equally spaced baffles. 
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(9) 


m 


./W 


k— tan k— L 
m m 


E C q A yqm 

w k tan k & 

y>q q q 


where 


yqm 


■/' 

yw 


(y+l)w qff(y Q " 


cos 


w 


cos 


m7Ty ( 

— 


dXr 


The integral for A,, 

yqm 

pactly as shown. 


was evaluated analytically but is written more com- 


The predicted frequencies obtained from numerical solution of Eq. 9, for 
the first- transverse mode in the main chamber, are shown in Fig. 1, along 
with some unpublished acoustic modeling data. These frequency results are 
not entirely unsatisfactory despite being somewhat low; however, the pre- 
dicted pressures were entirely unsatisfactory. Calculated pressures on 
each side of the interface are shown in Fig. 2. Because of the large dis- 
parity shown in this figure, the analysis and calculations were carefully 
checked for errors but none were found. Thus, it was concluded that a 
better approximation was needed. 


BAFFLE-COMPARTMENT APPROXIMATION 


The second approximation was developed from an assumed pressure profile of 
cos mrry/W on the compartment side of the interface. A somewhat more com- 
plicated characteristic equation was obtained (see Appendix B) : 


| —2. (k tan k 1) A - A 

v !e < M “ q q ^ M 

i-J W k tan k L 

m mm 


E 


£ 

—3- (k tan k &) 
w q q J 



0 


( 10 ) 
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NONDIMENS I ONAL FREQUENCY, wW/c 


MAIN CHAMBER APPROXIMATION 
FIRST TRANSVERSE MODE 


TWO-DIMENSIONAL CHAMBER 
WITH ONE OR FOUR BAFFLES 



0 0.1 0.2 0.3 0.4 0.5 


FRACTIONAL BAFFLE LENGTH, 1/ {l+l) 

Figure 1. Predicted Frequency Dependence on Baffle Length 
From the Main Chamber Approximation 
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PRESSURE (ARBITRARY SCALE) 


MAIN CHAMBER APPROXIMATION 
FIRST TRANSVERSE MODE 



Figure 2. Predicted Pressure Profiles on Each Side of the Main Chamber- 
Compartment Interface From the Main-Chamber Approximation 
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The pressure distribution on the compartment side of the interface for the 
til 

y compartment was given by 


- E^cos an 

q 

which is merely a Fourier series representation of cos muy/W. When numerical 
solution of Eq. 10 was attempted, no solutions were found to exist unless 
only the first term in q was retained, in which case all coefficients in 
the Fourier series (Eq. 11) are zero except the first. Thus, the pressure 
match at the interface was poor in this case as well; exemplary results are 
shown in Fig. 3. However, the frequency results improved. 

The iterative calculations were made with this compartment approximation 
as a starting condition, i.e., a zeroth iteration. 

ITERATIVE APPROXIMATION 


The characteristic equation obtained in this case (Appendix B) is the same 

as that shown in Eq. 10 except the coefficients A _ are replaced by an 
_(n+l) Uqm 


yq 


where 


and 


t (n+l) 

yq 


m 


m 

W 


y»q 


£ 

w 


(k tan k £) 

q q 


A a (n:)| 
yqm yq | 


k tan k L 
m m 


yqm 


( 12 ) 


,( 0 ) 

yq 


= A 


yqm 


The frequencies obtained by numerical solution of this characteristic equa- 
tion for the first- transverse, main-chamber mode are shown in Fig. 4, along 
with the acoustic modeling data. The curve denoted by "zero iterations" is 
the result from the compartment approximation just described. Agreement 
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PRESSURE (ARBITRARY SCALE) 


COMPARTMENT APPROXIMATION 
FIRST TRANSVERSE MODE 


TWO-DIMENSIONAL CHAMBER 


WITH ONE BAFFLE 



0.2 0.4 0.6 0.8 1.0 


TRANSVERSE POSITION, y/w 

Figure 3. Predicted Pressure Profiles on Ea’ch Side of the Interface 
From the Compartment Approximation 
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NOND I MENS IONAL FREQUENCY, wW/c 



FRACTIONAL BAFFLE LENGTH, £/U+L) 


Figure 4, Predicted Frequency Dependence on Baffle Length From the 
Iterative Approximation 
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between predicted and measured frequencies is quite good. The data show 
a somewhat higher trend at the longer baffle lengths, but this is believed 
to be due to inter comp artmental leakage in the acoustic model. 

The modeling data were obtained from a two-dimensional rectangular model 
(20 x 30 centimeters) with one or four baffles of (approximately) equal 
length and spacing mounted in one end. The largest share of these data 
(those obtained with four baffles) were obtained with relatively loose- 
fitting baffles. This lack of snug fit is likely to have caused some 
error in experimental results, especially with long baffles, because of 
intercompartmental • leakage. A frequency depression similar to that shown 
by these data has been observed repeated in engine tests and in other 
modeling tests (for example see Ref. 5) . 

The' characteristic equation as written in Eq. 10 and previously in Eq. 4, 
represents an overall energy balance. However, it appears more appropriate 
to consider an individual balance for each compartment. If this is done, 
the characteristic equation is similar to Eq. 10 except one summation over 
y is removed from each term. Fortunately, calculations made in this manner 
yielded the same frequencies as the overall equation. Thus, it was not 
necessary to worry about selecting the best "average" frequency. 

A satisfactory pressure match across the interface was obtained from the 
iterated calculation. Typical results are shown in Fig. 5 and 6 for 
chambers containing one and four baffles, respectively. Note that the 
error appears to diminish as the number of baffles contained in the chamber 
is increased. Hopefully, with a large enough number of terms in each 
series, and with sufficient iterations, exact agreement in the pressures 
can be obtained. 

A limited attempt to show a trend toward such convergence was somewhat 
discouraging, however. The average pressure difference (rms) between 
the two sides of the interface was calculated for one baffle configuration. 
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PRESSURE (ARBITRARY SCALE) 


FIRST TRANSVERSE MODE 


TWO-DIMENSIONAL CHAMBER 



0 0.2 0.4 0.6 0.8 1.0 


TRANSVERSE POSITION, y/W 

Figure 5. Predicted Pressure Profiles for a Single Baffle From the 
Iterative Approximation 
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PRESSURE (ARBITRARY SCALE) 



Figure 6. Predicted Pressure Profiles for Four Baffles From the 
Iterative Approximation 
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while the number of terms in m and q and the number of iterations were in- 
dividually varied from a base combination. The configuration was a single 
3- centimeter- long baffle in one end of a 20- by 30- centimeter chamber. The 
results (Fig. 7, 8 and 9) show that the errors are reasonably low, but 

they exhibit a disturbing trend toward increased error with increased q 
and, to some extent, with the number of iterations. However, these effects 
have not been properly investigated; convergence may indeed occur with a 
simultaneous increase in each of these terms. Further, the results may be 
affected by the small number of terms retained, which was necessary because 
of the small capacity of the time-sharing computer used to solve the equation. 

Based on the iterative method, pressure profiles were calculated at several 
longitudinal positions in the chamber. The results are shown in Fig. 10 
along with acoustic modeling data. The calculated profiles are symmetrical, 
whereas the measured profiles are not, probably because the baffle was not 
exactly centered and the model was driven from one side (x = 1.5 cm and 
y = 0) . Nevertheless, the agreement is regarded as good. 

Little attention was given to the possibility of calculating baffle com- 
partment modes. However, such a mode was indicated in some of the calcu- 
lations. The results are not presented here because the pressure match 
at the interface was poor. 

TWO-DIMENSIONAL CHAMBERS WITH GAIN/LOSS 
BOUNDARY CONDITIONS 

The effects of pressure- coup led gains and losses were simulated to a degree 
by adding nonzero acoustic admittance values at both ends of the chamber. 

A loss- type condition was imposed at the nozzle end and a gain- type condi- 
tion at the injector end. These effects were added so that the influence 
of baffles on (pressure-coupled) stability could be investigated. Because 
no steady flow was included in the calculation, effective or virtual boundary 
conditions were used to compensate for the omission. However, this distinc- 
tion between actual and virtual boundary conditions only affects the choice 
of numerical values to be examined and, hence, is of little importance. 
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RMS ERROR, Ap/p RMS ERROR, Ap/p 

max max 


FIRST TRANSVERSE MODE 



0 2 4 6 8 10 12 

MAXIMUM TERM RETAINED IN SUMMATION OVER m, i.e., m 

Figure 8. Dependence of Error on Number of Terms in m Retained 



02 46 8 10 12 

MAXIMUM TERM RETAINED IN SUMMATION OVER q , i.e., q 

Figure 9. Dependence of Error on Number of Terms in q Retained 
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ITERATIVE APPROXIMATION 
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The boundary condition was written in terms of a specific acoustic 
admittance: 


y = 



( 13 ) 


To partially compensate for steady through flow at the boundaries, this 
admittance may be replaced by (Ref. 6 and 7) 

pc^T = YM je[ (14) 

where M is the steady-flow Mach number through the boundary and y/£ is the 
response function (fractional perturbation in mass flux divided by the 
fractional perturbation in pressure) . Equation 14 was used only to estimate 
numerical values, for the acoustic admittances. 


The iterative form of the characteristic equation may be rewritten to include 
these boundary conditions (see Appendix B) . The result is 


m 


m 

W 




—S- k tan(k Z + ip ) a^ n+ ^ A I 
w q q r q yq yqml 


k tan (k L + l p ) 
m m m 


y>y 


-3- ia^ n+ ^l k tan(k Z + ip ) = 0 

» yq q q q 


(15) 


where 


tan = -J^y N /k m 


tan ib = -jky T /k 

q J q 
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and y^j and y^ are the specific admittance values at the nozzle and injector 
ends, respectively. 

The effects of baffles on the stability of this chamber were investigated 
by comparing the stability requirements with and without a baffle. For a 
fixed nozzle loss, the injector-end admittance necessary for neutral stabil- 
ity was calculated, thus defining a stability limit curve. This was done 
by solving Eq. 15, with a numerical root-finding technique, for the y^ that 
corresponded to particular, real eigenvalues and a constant nozzle admittance. 


A stability limit curve was also calculated for the unbaffled case from 


j -F ta "< k m L * V - Vi 


Ci«) 


The nozzle admittance was estimated from that for a zero-length nozzle with an 
entrance steady-flow Mach number of 0.266 and y = 1.25. The response func- 
tion for a zero- length nozzle is (y + l)/2y; thus, a value of y^ = 0.3 was 
chosen. 

Stability limit curves for the unbaffled case and for a single baffle are 
shown in Fig. 11. As defined, the admittance has a driving effect if the 
real part of the admittance .is negative. Therefore, these calculations 
indicate the engine is less stable with a baffle than without. 

This result is not very surprising in view of the observed effect of baffles 
on the pressure distribution within the chamber; the baffles increase the 
amplitude at the injector end relative to the nozzle end. According to 
Cantrell (Ref. 6), neutral stability is obtained when 

f\ p 2 Re | y| dS = 0 (17) 

If the amplitude at the injector end is increased, the real part of the 
admittance at that end must be reduced to balance the same nozzle loss. 
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STABILITY LIMIT CURVES 
FIRST TRANSVERSE MODE 
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This implies that baffles improve the stability of an engine in some other 
manner than reducing pressure coupling, at least to the extent that pressure 
coupling is approximated by an acoustic admittance. 
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CONCLUSIONS 


Results from this study indicate that a largely satisfactory mathematical 
technique has been developed for calculating the wave motion in baffled 
chambers. This technique has been applied to several simple cases with 
good success, but additional work is needed to develop a thorough under- 
standing of the oscillatory behavior of baffled chambers. The ability to 
obtain an essentially continuous pressure profile along the baffle tips has 
been demonstrated with this technique. Requirements for continuity of 
velocity and energy flux at the interface between the baffle compartments 
and main chamber are automatically satisfied. 

The effects of a single baffle on the stability of a chamber with nonrigid 
walls (gain/loss- type boundary conditions) has been successfully analyzed 
for one particular geometry. Thus, the ability to generalize the method 
for nonzero boundary conditions has also been demonstrated. Convergence 
of the iterational scheme has not been completely demonstrated. However, 
very good numerical results were obtained and it is likely that adequate 
convergence can also be obtained. Convergence appears even less restrictive 
when the chamber contains more than one baffle. 
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APPENDIX A 


SECULAR EQUATION TECHNIQUE 

Morse* suggests approximate solution of integral equations, such as those 
encountered herein, by assuming an eigenfunction expansion for the dependent 
variable. A secular determinant is thereby obtained. This approach is more 
appealing in some respects than the variation-iteration technique described 
previously because no initial estimate of the solution is required. For 
this reason, considerable emphasis was placed on this technique during the 
program. However, satisfactory numerical results were never obtained. The 
reason for this failure is not known. The method is believed applicable; 
therefore, the failure may be due to a mathematical or programming error, 
although both the analysis and the computer program have been checked 
thoroughly. 


Because the secular equation method is still believed applicable, it will 
be briefly outlined here. Both the Green's functions and the normal gradient 
(?) can be expanded in terms of a convenient set of eigenfunctions; the 
Green's functions can be written explicitly, but the coefficients in the 
expansion for ? are left unspecified. The Green's functions appear as 



v' vV 

Z - J 2 2 

n VA n (n n -k 2 ) 


and ? as 


5t? <P " £ c m Vfy 

m 


(A-l) 


(A- 2) 


*Morse, P. M. , and K. U. Ingard: Theoretical Acoustics , McGraw-Hill Book 

Co., Inc., New York, 1968. 

Morse, P. M. , and H. Feshbach: Methods of Theoretical Physics , McGraw- 

Hill Book Co., Inc., New York, 1953. 
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Insertion of Eq. A-2 into the set of integral equations, Eq. 4, and 
employing the orthogonality properties of the eigenfunctions, leads to 
a set of linear algebraic equations in c , each equation being equal to 
zero. Therefore, the determinant of the coefficients of the c array 
must equal zero to give a nontrivial solution. This determinant is a 
secular determinant; it represents a characteristic equation that speci- 
fies the eigenvalues (frequencies) of the combustion chamber. 

In principle, the secular determinant is of infinite order, but the indi- 
cations are that it can be truncated at some reasonable size. Once the 
eigenvalues are determined, the array of equations in c can be solved for 

each of the coefficients, each c , in terms of one of them. The function 
^ m 

£(Tq) is thereby specified and the pressure distribution throughout the 
combustion chamber may be calculated. 

As noted above, the secular- equation formulation was developed and programmed 
for numerical solution. Initial results were promising, but eventually it 
was found that unacceptable frequency and pressure results were being ob- 
tained. Specifically, a pressure discontinuity was predicted about midway 
between baffle tips. The reason for this was not found. 



APPENDIX B 


DEVELOPMENT OF EQUATIONS 

In this appendix the mathematical details concerning the variational solu- 
tion of the wave equation (actually, the Helmholtz equation) for closed, 
two-dimensional, baffled chambers are summarized. The initial section is 
concerned with development of the variational form of the characteristic 
equation, while the second section deals with the development of the ap- 
proximate velocity distributions necessary for the main-chamber and baffle- 
compartment approximations. The iteration scheme is also developed there- 
from, and finally, the equations for the nonzero admittance case are 
summarized. 

VARIATIONAL FORMULATION 

The variational formulation used during this program is patterned closely 
after that described by Morse and Ingard* . Employing the variational func- 
tion developed by them, with a slight generalization for multiple cavities, 
a characteristic equation was obtained. 

/«/ G 5 dS 0 ^ 5 iS 0 ^ - 0 ( B '« 
t t t p 

where E, is proportional to the normal component of velocity at the interface 
between the main chamber and the compartments. This equation was obtained 
from the variational method for an approximate normal velocity of the form 
u = A?, the value of the amplitude (A) being optimized by the method. 

Equation B-l may be rewritten as 

/(Pa - P bp > dS = 0 (B-2) 


*Morse, P. M. , and K. U. Ingard: Theoretical Acoustics , McGraw-Hill 

Book Co., Inc., New York, 1968, p 680. 
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because, except for an arbitrary amplitude coefficient. 


P a (?) = -jpio/'l G a (r|rjj) ?(r|j) dS Q (B-3) 

S t 

P by (?) = jpo3/'G by (?|^) ?(rg) dS 0 (B-4) 

Equation B-2 implies that the best estimates, in the variational sense, of 
the eigenvalues (frequencies) of the baffled chamber are obtained by requir- 
ing continuity of the energy flux at the interface. Thus, requirements for 
continuity of energy flux and velocity at the interface are automatically 
satisfied. However, calculated, pressures on each side of the interface will 
not be equal unless the exact normal velocity distribution is used. 

Clearly, the next task is to develop a means for finding suitable estimates 
of the normal velocity distribution. Rather than attempting to choose a 
velocity distribution directly, however, approximate pressure distributions 
were selected and the corresponding velocity distributions were calculated. 

MAIN-CHAMBER APPROXIMATION 


This approximation was obtained by assuming a pressure profile on the main 
chamber side of the form' 


P a (L,y) = cos — £ 

However, 

P a ( L >y) = A c?S|? o } ^ dy o 


(B— 5) 


(B-6) 
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and 


<y?l? 0 ) - 


-E 


miry mT1 ^0 

jn cos W cos — 

W k sin k L 
m m 


F m« 


where 


F m (x) = < 


cos k m (x Q - L) cos k ffl x x < x Q 


cos k x n cos k (x - L) x > x. 
m 0 m 0 
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Therefore, 
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3p a m7T y c 

77 COS 7T- 

3x o w 


dy r 


k tan k L 
m m 


W/e- 

m 


m = m 


m / m 
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Further, from the integral expression for p , 

cL 


ap a a,y) v ^ 
w 


Bx 


cos 


mny 

W 


fK 

J 3x n 


a mTTy 

ax" cos w % 
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Consequently, from Eq. B-8 and B-9, 


9p a ( L >y) 

dx 


= E, = -k_ tan k__L cos 

^ m m W 


(B- 10) 


Introducing this expression into the characteristic equation (Eq. B-l) 
leads to 


&-/W 

nr + 

k— tan k— L 
m m 



k 
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Pqm 


tan k l 

q 


0 


(B-ll) 
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where 


/ (y+l)w m7ry 0 qiT(y 0 - yw) 

cos — 7 -. — cos — dy„ 

W w 1 0 

■ 1 yw 


COMPARTMENT APPROXIMATION 

The compartment approximation was obtained in a similar manner by 
a compartment- side pressure profile of 

Pby (L ’>° = cos 


for 


yw < y < (y + 1) w 


th 

The Green's function for the y*" 1 compartment may be written as 




E COS 



w 


q"ty- W) ^ 0^0- 


W 


W 


k sin k i 

q q 


F q (x) 


where 


F q (x) * 


cos k (x„ - L,.) cos k (x - L) x < x„ 
q 0 t q J 0 


cos k (x„ - L) cos k (x - L.J x > x„ 
q 0 q t 0 


L = L + l 


For this approximation 


9p. 


_by 

9x 


Z — 3- (k tan k £) 
w q q 


, A _ cos 

q Uqm 


gn (y-yw) 
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assuming 


(B- 12) 


(B- 13) 
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Insertion of this approximation into the characteristic equation (Eq. B-l) 
leads to Eq. 10. The corresponding pressure expressions were obtained 
from Eq. B-3 and B-4. 

ITERATION SCHEME 

Neither of the last two approximations was considered adequate. Therefore, 
an iteration scheme was developed to improve an initial estimate. A pressure 
distribution was assumed of the form 


PbpCi-.y) - 


a („) cos t u fr .. - ,, uw) . 
yq 


w 


(B-15) 


for the interval yw < y < (y+l)w. 


The corresponding pressure gradient is 


8 V L,y) = V !a k tan k ^ cos Vir -.jy . ) 

Bx x — ' w pq q q w 

q 
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Further, the corresponding pressure expression for the main chamber side is 


p a ( L >y) = 


-z 


cos 


miry 

W 


m 

W k tan k L 
m m 


y,q 


a*" 11 "* (k tan k A 
w yq q q yqm 
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Now, by equating the right side of Eq. B-17 to 311 improved estimate 

may be obtained. From this, the following expression was obtained: 


(n+1) 

yq 
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(Z 

\y,q 


£ 

-5. a 
w yq 


(n) 


k tan k SL 

q q . 


k tan k L 
m m 


yqm 
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The characteristic equation obtained in this case is similar to Eq. 10, 

the only difference being a replacement of A —in that equation by the 
(n+1) , qm 

a shown above. 

uq 

NONZERO SURFACE ADMITTANCE 

The influence of nonzero acoustic admittance boundary conditions at the 
closed ends of the chamber may be readily included in the analysis by 
choosing Green's functions which satisfy the same boundary conditions as 
those satisfied by the pressure. The appropriate Green's functions are 
then 


where 


and 


G a t? l ? o ) = £ 


221 ... my o 


COS — rf- COS 

m W W 


_ W k sin(k L + ip ) m 
m m "• m r nr 


F fx) 
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F (x) = 

m J 


cos k (L - x.) cos(k x + ip ) 
m v Cr ^ m r m 


cos (k x« + 4> ) cos k (L - x) 
m 0 nr m^ 


x < x. 


x > x r 


tan *m * 
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cos 


qiT(y- w) ^ &0 ^ 

- 1 — cos 

w w 

k sin(k H + ip ) 

q q q 


F q (x) 
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where 


F q (x) 


cos [kq(L t - x q) + cos k^(x - L) x < x 0 
cos k^(x 0 - L) cos[k q (L t - x) + Jp J ] x > x Q 
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and 


tan 


s k V k q 


When these Green's functions were employed, expressions for pressure and 
the characteristic equation were obtained which are very similar to the 
iterative expressions just described, except that the terms and k^L 
are replaced by k^& + iJj and k^L + ip^, respectively. That is, the char- 
acteristic equation is 
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u,q 


V"* -3- k tan(k Z + ip ) a^ n+ ^ A 
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Uq 


Vqm 


k' tan (k L +" \p )' 
m mm 
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k tan(k Z + ip ) 

q q V 


(B- 21) 
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APPENDIX C 


COMPUTER PROGRAMS 

This appendix summarizes the four computer programs used during this study 
Philosphy and derivation of the working equations were described previously 
in the main body of this report. This appendix describes the numerical 
methods used to solve these equations, and includes input, output, and opera- 
tional information as well. The nomenclature of all four programs is closely 
related and is presented in the Nomenclature subsection. 

The four computer programs used in this study are: 

1. SSMCE4--solves the characteristic equation for the rigid boundary 

case 

2. SSBPD2--calculates the pressure on the main chamber side of the 

interface across the baffle tips 

3. SSVEL2-- calculates the pressure on the compartment side of the 

interface across the baffle tips 

4. SSMCE8--solves the characteristic equation for the active boundary 

case 

These four programs were run on a timesharing computer system (General 
Electric 420 computer). All four computer programs are not written effi- 
ciently with respect to running time; some statements could be eliminated 
without inhibiting the calculations. This is especially true for programs 
SSBPD2, SSVEL2, and SSMCE8. The reason for the inefficiency is that the 
main portion of each of these was pirated from program SSMCE4. 

SSMCE4 

Computer program SSMCE4 was written to solve the iterative characteristic 
equation for the rigid boundary case. The method used to solve this 
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equation was to calculate the value of the function while incrementally 
changing"phee" (cf>) until the value of the function changed signs. The usual 
practice was to calculate the characteristic equation over a large interval 
by taking large increments of "phee." After noting the interval in which 
the value of the function changed signs, that interval was subdivided into 
smaller increments. The procedure was repeated until an accurate value of 
PHEE was obtained. 

The four subroutines used to calculate repetitive terms are: 

1. Subroutine TAN(a,b,c) = k. tan k.£ 

2. Subroutine TANH(a,b,c) = -k^ tanh k^ib 

3. ARG(a,b,c, x) = A 

4. ARI (a,b,c, x) = 

The computer program names for the various parameters in the characteristic 
equation are presented in the Nomenclature subsection of this appendix. 

The input information required is: 


1 . 

PHEE1 

= 

Starting phee value 

2. 

ICON 

= 

Index for transfer 

3. 

NUMB 

= 

Number of points (phee values) at which the function 
will be evaluated 

4. 

STEP 

- 

Increment for increasing phee from the initial vajme 

5. 

MU1 

= 

Index of the baffle compartment being examined 

6. 

NS 

= 

Number of iterations on a, 

pq 

7. 

Q1MAX 

= 

Number of q+1 terms taken in the a iteration 

H yq 

8. 

MUMAX 


Total number of baffle compartments 

9. 

QUMAX 


Total number of q+1 terms taken in the evaluation 
of the characteristic equation 
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10. 

MMAX 

= Total number of m+1 terms taken in the 
of the characteristic equation 

evaluation 

11. 

MBAR 

= Index of the particular acoustic mode 

to be investigated 

12. 

SMLW 

= Width of a baffle compartment 


13. 

WDTH 

= Width of the main chamber 


14. 

SMLEL 

= Length of the baffles 


15. 

XLNGTH 

= Total length of the chamber 



Output from SSMCE4 includes the phee value and the corresponding value of 
the characteristic equation. 

SSBPD2 

This program calculates the pressure across the baffle tips on the main 
chamber side of the interface, and was written for Eq.B-17. With the phee 
value fixed for a given set of parameters (solution from SSMCE4) , Eq. B-17 
is an explicit relationship between pressure and the position coordinate 
(y) . Pressure at the baffle tips was calculated across the width of the 
chamber. 

Input for SSBPD2 is the same as that for SSMCE4. The phee value input for 
SSBPD2 is the root obtained from SSMCE4 for a given set of parameters. 
Output from SSBPD2 includes pressure and the corresponding values of y. 

SSVEL2 


This program was used to calculate the pressure at the baffle tips on the 
compartment side of the interface. Equation B-15 describes this pressure 
Again, for a given phee value, the pressure is a function of the position 
coordinate (y) . Thus, for a fixed set of parameters, the pressure can be 
readily calculated for a given y. 
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Input to SSVEL2 is the same as that for SSBPD2, with the following additions: 


y = Initial starting point of y 

Step = The incremental step in y 
SMLEL1 = SMLEL (nonfunctional parameter) 

Output includes pressures (p^ ) and the corresponding position coordinate (y) . 

SSMCE8 

This program was written to determine the stability limit of a combustion 
chamber with active boundaries located at both ends. The characteristic 
equation is in complex notation, and the root of the characteristic equation 
is complex eigenvalue (the real part is the nondimens ional frequency and 
the imaginary part is the nondimens ional damping coefficient) . The charac- 
teristic equation for the active boundary case is given by Eq. B-21. Thus, 
for a given set of parameters, nozzle admittance, and injector admittance, 
the complex root which satisfies the characteristic equation specifies the 
frequency and the damping coefficient of the system. However, for the 
stability limit of the system, another route can be taken. Because the 
damping coefficient is zero, the above equation can be solved for y rather 

than <j>, with fixed values of y and phee (REAL) . The acoustic admittance 

N 

thus obtained defines the maximum amount of acoustic energy (related to the 
admittance) that can be pumped into the system and still have the system 
stable. This avenue was taken to evaluate the stability limit. 

The root-finding technique used to determine y^ was Newton's method, with 
a finite-difference approximation for the derivative. 

With this approximation, a computer program was generated from SSMCE4 with 
a few alterations. Convergence was not rapid but was sufficient. 
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The program was written using complex algebra. The input parameters for 
SSMCE8 included the input for SSMCE4 and admittance. However, y , y^., and 
phee are complex numbers in this program. Output includes the y^'s and 
the corresponding values of the characteristic equations. 
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NOMENCLATURE 



SSMCE4 

SSBPD2 

SSVEL2 

SSMCE8 1 

£ 

m 

EPM 


EPM 


EPM 


EPM 


m 

M 


M 


M 


M 


w 

WDTH 


WDTH 


WDTH 


WDTH 


U 

MU 


MU 


MU 


MU 


q 

QU 


QU 


QU 


QU 



EPQ 


EPQ 


EPQ 


EPQ 


w 

SMLW 


SMLW 


SMLW 


SMLW 


k tan k £ 

q q 

Subroutine 

TAN 

Subroutine 

TAN 

Subroutine 

TAN 

Subroutine 

TAN 

a qq 

Subroutine 

ARI 

Subroutine 

ARI 

Subroutine 

ARI 

Subroutine 

ARI 

A 

yqm 

Subroutine 

ARG 

Subroutine 

ARG 

Subroutine 

ARG 

Subroutine 

ARG 

L 

XLNGTH 


XLNGTH 


XLNGTH 


XLNGTH 


£ 

SMLEL 


SMLEL 


SMLEL 


SMLEL 


<P 

PHEE1 


PHEE1 


PHEE1 


*PHEE1 


7T 

It 

PHEE 


PHEE 


PHEE 


*PHEE 


TT 

PI 


PI 


PI 


PI 


Pa 



FIRSUM 






. y 



^coordinate 

y- coordinate 



Pby 





QUSUM 




y q 







*YQ 


y m 







*YMA 



* Comp lex numbers 
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SSMCE4 


09J07 R0CKET 101969 


' KEY 

50C CHARACTERISTIC EQ ' N F0R THE BAFFLE PROBLEM 

55 INTEGER QU>QQjQUMAXjQU2>Q1MAX 

56 DIMENSI0N RET9C20.»20.*1 )iXl <20-»20 j1 ) 

57 

60 10 PRINT 915 

65 91 5 F0RMATC1HO*2X*26H1NPUT; PHEE* ICON* NUMB* STEP) 

70 INPUT jPHEEI * ICON* NUMB* STEP 

71 PRINT 917 

72 917 F0RMAT<1HO>2X>19HINPUT: Mljl * NS*Q1 MAX > 

73 I NPUT * MU 1 » NS* Q 1 MAX 

75 IF ( IC0N12O* 20* 30 

76 30 PRINT 916 

77 916 F 0RMAT C 1 HO * 2 X* 51 HINPtJT :MUMAX»Q'JMAX*MMAX* MBAR* WDTH* 
78&SMLW*SMLEL*XLNGTH > 

7 9 INPUT*MljMAX*QUMAX*MMAX*MBAR*WDTH*SMLW* smlel*xlngth 

80 PHEE=PHEE1 /WDTH 

81 20 PRINT 910 

82 N= 1 

83 PRINT 91 1 

84 PRINT 918*PHEE1*WDTH*SMLW*SMLEL*XLNGTH 

85 PRINT 912 

86 PRINT 91 9*MMAX*MUMAX*QUMAX*MBAR 

87 PRINT 913 

88 FMBAR=3 

89 15 00 170 M(J2= 1 * M(JMAX 

90 D0 171 QU2=1>QU«AX 

92 FMBAR=FL0AT(MBAR) 

93 Pl=3. 141592 

94 RET9(QU2*MU2*MBAR>=0 • 

95 171 CONTINUE 

96 170 CONTINUE 

97 00 172 MU2= 1 * Ml JMAX 

98 Q.U2=1 

99 FMU2 = FL0ATCMy2-*l ) 

101 FQU2=FL0AT(QU2-1 ) 

102 IF ( F MBAR )21 *21 >22 

103 22 IFCFQU2>31*31*32 
10421 IF (F QU2 1 33* 33* 34 

105 33 RET9<QU2*Mfj2*M6AR)=SMLW 

107 60 T 0 11 

108 34 RET9<QU2*MU2*MBAR>=0 -0 

109 G0 T 0 1 1 

110 PRINT 911 

111 31 XI=FM8AR*PI/WDTH 

1 12 RET9CQU2*MU2*MBAR)=CSINIXI*CFMU2+1 .0)*SmLW) 

1 13&-SINCXI*FMU2*SMLW))/XI 

114 G0 T 0 11 
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115 32 JK=FQU2 

1 16 XA=(-1 .0)**JX 

117 XF=FQU2*Pl/S«LW 

118 xb=fmbar*pi*$mlw/wdth 

119 XC=XB*FMU2 

120 PRINT 912 

121 XU = XB*CFMij2+l -0) 

122 XE=FMSAR*PI/WDTH 

123 IFCXE-XF)58, 59, 58 

124 59 RET9CQU2,M!j2,MBAR)=SMLW*C0SCFMU2*FQlj2*Pl ) /2»0 

125 G0 T 0 11 

126 58 XG=XE/<XF*XF-XE*XE) 

128 RET9CQU2, MU2, MBAR ) =XG* C S I N C XC ) - XA*SI N< XO > > 

129 11 CONTINUE 

131 173 CONTINUE 

132 172 CONTINUE 

133 CALL ARI<PHEE,WDTH,SMLEL, SMLW,xLnGTH,RET9, 

1 34&MUMAX,QUMAX,MMAX,M8AR, Mijl , NS, Q1 MAX, FMljO, F QUD, X 1 > 

135 Pl=3. 141592 

140 FIRSIJM=0.0 

145 D0 200 M=1,MMAX 

150 FM=FL0AT<M-1 ) 

155 AA=FM*PI/WDTH 

160 DISCR=PHEE*pHE£-AA*AA 

165 IF C D I S CR ) 50 , 51 >52 

170 50 ART=SQRT(-DISCR)*XLN6TH 

175 I F CART - 4 • ) 8 50 ,851,851 

180 850 CALL TANHY < D I SCR* XLNGTH , RET > 

185 G 0 T 0 852 

190 851 RET=-SQRTCABSCOISCR) ) 

195 852 60 TO 300 
200 51 PRINT 980 
205 G0 TO 10 

210 52 CALL TANCOISCR, XLNGTH, RET) 

211 300 3MUSUM=0.0 

212 D0 220 MU = l»MiJMAX 

213 FMU=FL0ATCMU-1 » 

215 QSUM=0.0 

216 CSUM=0*0 

220 DO 210 QU=1,QUMAX 
225 FQIJ = FL0ATCQU-1 ) 

230 IF CFQU >60,60,62 

235 60 Ep Q U = 1 *0 

240 G0 T 0 301 

245 62 EP0U=2.0 

250 301 AB=FQU*PI/SMLW 

255 DISCr2=PHE£*PHEE-AB*AB 

260 IF CDISCR2 >65,66,67 

265 65 ARP=SQRT(-0ISCR2 )*SMLEL 

270 IFCARP-4. >860,861,861 

275 860 CALL TANHY CD I SCR2 , SML EL, RET2 > 
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280 G0 T 0 347 

285 861 RET2=-SQRTCA8SCDISCR2)> 

286 G0 T0 347 
290 66 PRINT 980 
295 G0 T0 10 

300 67 CALL TANC0ISCR2* SMLEL*RET2 ) 

330 347 FM8AR=FL0AT CMBAR ) 

331 CALL ARGCFQU*FMU*FM*SMLW* WDTH*RET3> 

333 FQR=0 »0 

336 FMU0«FMU 

337 FQUD=FQU 

341 QSUM=QSUM+EpQU*RET2*RET3*Xl CQU* t v iU* M8AR ) /SML vi 

342 

344 FMU1=FL0AT(MU1 -1 ) 

345 210 CONTINUE 

346 IFCFMU-FMU1 )4* 5* 4 

347 5 DSUM=QSUM 

355 4 SM!JSUM=SMUSUM+QSUM 

360 220 CONTINUE 

365 IF CFM )70*70*71 

370 70 EPM=1.0 

375 G 0 T0 303 

380 71 EPM=2.0 

385 303 FIRSUN=FIRSUM+SMUSUM*DSUN*EPM/(RET*WDTH) 
390 200 CONTINUE 
400 SECSUM=0»0 

405 MU=MU1 

406 FMU=FL0ATCMU) 

410 FMU=FL0ATCMU-1) 

415 FISU0U=0.0 

420 00 240 QU=l*QUMAX 
425 FQU=FL0AT<QU-1 ) 

430 IF (FQ U ) 80* 80* 8 1 

435 80 EPQU=1 »0 

440 60 T0 500 

445 81 EpQU-2 »0 

450 500 AC=FQU*PI/SMLW 

455 DISCR3=PHEE*PHEE-AC*AC 

460 IFCDISCR3)90*91 *92 

465 90 ARS=SQRT<-DISCR3>*SMLEL 

470 IF < ARS"4» ) 400 * 40 1*401 

475 400 CALL TANHY CD I SCR 3* SML EL* RET 5) 

480 G0 T 0 410 

485 401 RET5=-SQRTCA8S<DISCR3> > 

486 G0 T041O 
490 91 PRINT 980 
495 G0 T0 10 

500 92 CALL TAN(DISCR3*SMLEL*RET5) 

579 410 FQp=0-0 

583 FM(JD=FMU 

584 FQUO=FQU 
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600 FISUQU=T ISUQU+KET5*X1 <aU*MU*M8AR)*Xl < QU* «U* MbAR ) * 

601&EPQU/SMLW 

605 240 CONTINUE 

615 SECSUM=SECSUM+F ISUQU 

617 953 F 0RMAT <F 12 • 6) 

626 T0SUM=FIRSUM+SECS(JM 

627 PHEE1 =PH£E*WDTH 

630 PRINT 952*N*pHEEl *T0SUM 

635 IF < N- NUMB )975*975*976 

640 975 N=N+1 

645 PHEE=PHEE+STEP 

650 60 T0 15 

655 976 60 T0 10 

660 910 FORMAT (35HINITIAL VALijES F 0R THIS CALCULATION//) 

665 911 F0RMAT<1H0*3X*4HphEE*4X*4HWDTH*4X*5HSMLW *3X*5HSMLEL* 
666*2X*6HXLN6TH) 

670 918 F0RMATCF8»4#4F8*2//> 

675 912 F0RMATUHO*4X*4HMMAX*3X*5HMUMAX*3X* 5HQIJMAX* 4X* 4HM8AR ) 
680 919 FORMAT C4C5X*I3)///) 

£81 935 FORMAT ( 1 HQ*2X* 6HQQSUM=* F 1 2 • 5* 2X* 3HQQ = » 12* 2 X* 3HQiJ= * 
682&I2*2X* 3HMU=* 12) 

683 936 F0RMATClHO*2X*7HF > lSUaU=*F 12.5*2X*3HQa = * I2*2X*3HQU=* 
684&I2*2X*3HMU=* 12) 

685 913 FORMAT ( 1 HO* 1 X* 1 HI * 5X* 4HPHEE* 6X* 6HFEF UNC ) 

686 937 F0RMAT<1HO*2X* 7HSECSUM=* F 12 . 5>2X* 3HQQ=* I2*2X*3HQU=* 
687&I2*2X*3HMU=* 12 ) 

690 980 FORMAT C9HUISCR=0.0) 

693 932 FORMAT ( 1 HO* 2X* 7HSMUSUM= * F 12 . 5*2X* 3HMy=* I2*2X* 3HQ|J=* 
694&I2*2X*2HM=* 12 ) 

695 933 F0RMATC1HO*2X*5HQSUM = *F12.5*2X*3HM|J = * I2*2X*3HQU=* 
696&I2*2X*2HM=* 12) 

697 934 F 0RMAT ( 1 H0*2X* 7 HF IRSUM=* F 1 2 • 5* 2X* 3HMU=* 12* 2X* 3HQ(J = * 
698&I2*2X*2HM=* 12) 

700 952 F0RMATCI3>2<1PE12.4)) 

705 END 

710 SUBROUTINE TANHyCX*Y*Z) 

715 RR=SQRT CABS <X ) ) 

720 RS=RR*Y 

725 Z=-RR*CEXPCRS)-EXPC-RS) )/CEXpCRS)+EXpC-RS) ) 

730 RETURN 
735 END 

740 SUBROUTINE TAN<X*Y*Z) 

741 RT=SQRT<X) 

742 R|J=RT*Y 

743 Z=RT*CSINCRU)/C0S<RU)) 

745 RETURN 

746 END 

747 SUBROUTINE AR 6 <TA* TB* TC* TD* TE* TF ) 

748 Pl=3. 141592 

7 49 IF CTC ) 720* 720*730 
751 730 IFCTA>731*731*732 
7 52 720 I F ( TA) 7 33* 7 33* 7 34 
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753 733 TF = TD 

754 G0 T 0 710 

755 734 T F=0 • 

757 G0 T0 710 

758 731 XI=TC*PI/TE 

759 TF=CSIN«XI*CT8+1 .03+TD3 -SlN<Xl*TB*TD3 3 /XI 

760 G0 T0 710 

761 732 JK=TA 

762 XA=(-1 .03**JK 

763 XB=TC*Pl*TD/TE 

764 XC=XG*TB 

765 X0=XB*CT8+1 .03 

766 XE=TC*PI/TE 

767 XF=TA*pI/ro 

768 IF<X£-XF3760*761*760 

770 761 TF=TD*C0S<TA*T8*PI3/2.O 

772 G0 T0 710' 

773 760 XG=XE/CXF*XF-XE*XE3 

774 TF=XG*CSINCXC3-XA*SIM(X03 3 

776 710 CONTINUE 

777 RETIJKN 

778 END 

779 SUBROUTINE ARI CpHEE* WDTH* SMLEL* SMLW* XLNbTH* RET9* 
780AMUMAX*QUMAX*MMAX*MBAR*MiJl *NS*Q1 MAX* FMtjD* FQ(jB* X 1 3 

781 INTEGER QU*QQ*QUMAX*QU2*Q1MAX 

782 DIMENSION RET 9<20 *20* 1 3 * x 1 <20 *20* 1 3 

783 NN= 1 

784 1 DO 175 MU2=1*MUMAX 

785 DO 176 QU2=1*QUMAX 

786 Pl=3. 141592 

787 FIRSUM=0. 

788 DO 200 M= 1 * MMAX 

790 FM=FL0AT<M-1 ) 

791 AA=FM*PI/WDTH 

792 DISCR=PHEE*PHEE-AA*AA 

793 I F < D I SCR 3 50 * 5 1 * 52 

794 50 ART=SGRT<-DISCR)*XLNGTH 

796 IFCART-4. 03850*851 *851 

797 850 RR=SQRT CABS <D I SCR 3 3 

798 RS=RR*XLNGTH 

799 RET = -RR*<EXP<KS3-EXP<-RS> > / < EXP <RS> +EXP < -RS > > 

801 GO TO 852 

802 851 RET=-SQRT<ABS<DISCR>> 

803 8 52 GO T0300 

804 51 ST0P 

806 GO TO 10 

807 52 RT=SQRT<DISCR> 

808 RU=RT ♦XLNGTH 

809 RET=RT^<SIN<RU3/C0S(RU3 3 
611 300 SM|JSUM=0. 

812 DO 220 MU*1*MUMAX 

813 FMU*FL0AT<MU-1 3 

814 QSUM*0. 

815 D0 210 QU=1*Q1MAX 
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816 FQtJ=FL0ATCtfU) 

817 FQU=FQU-1 »0 

818 IF(FQU>60*60*62 

819 60 EPQU=1 .0 

821 60 T 0 301 

822 62 EPQ U=2 • 

823 301 AB=F0U*PI/SMLW 

824 DISCR2=PHEE*PHE£-AB*AB 

826 IF (D1 SCR2 >65*66*67 

827 65 ARP=SQRTC-DISCR2)*SMLEL 

828 IF ( ARP-4 »0 > 860*861*861 

829 860 RR = SQHT CABS CD I SCR2 ) ) 

830 RS=RR*SMLEL 

832 RET2=“RR*CEXP(RS)-EXP<-RS> )/(EXP(RS)+EXPC-RS) ) 

833 60 T0 347 

834 861 RET2=-S0RT<ABS<DISCR2>> 

836 60 T0 347 

837 66 PRINT 980 

838 10 STOP 

839 67 RT=S0RT(0ISCR2) 

841 RU=RT*SMLEL 

842 RET2=RT*CSIN<RU)/C0S(KU> > 

843 347 IF<FM>720*720*730 

844 730 IF CF0U>731*731*732 

846 720 IF(FQU>733*733*734 

847 733 RET3=SMLW 

848 60 T 0 710 

849 734 RET3=0 • 

851 60 T0 710 

852 731 XI=FM*PI/WDTH 

853 RET3=CSIN(XI*CFMU + 1 »0 > *SML W > -S I N< X I*F MU*SML W) ) /XI 

854 60 T0 710 

856 732 JK=FQU 

857 XA=(-1 .0)**JK 

858 XBsFM*PI*SMLW/WDTH 

859 XC=XB*FMU 

861 XD=XB*(FMJ+1 . ) 

862 XE=FM*pI/WDTH 

863 XF=FQU*PI/$MLW 

864 IF(XE-XF >760* 761*760 

866 761 RET3=SMLW*C0SCFMU*FQU*PI)/2.O 

867 60 T0 710 

868 760 X6=XE/CXF*XF-XE*XE> 

869 RET3=X6*CSIN(XC)-XA*SIN(X0) ) 

870 710 CONTINUE 

871 FQU£)=FL0AT<QU2-1 > 

872 FMUD=F L0AT <MU2- 1 > 

873 FMBAR=FL0ATCMBAR) 

874 QSUM=QSUM+EPQU*RET2*RET3*RET9CQU* MU^MBAR)/SNLW 

876 210 CONTINUE 

877 SMUSUM*SMUSUM+QSUM 

878 220 CONTINUE 

879 IF(FM)70*70*71 
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881 ?0 EPM= 1 • 

882 60 T0 303 

883 71 EpM=2.0 

884 303 IF<FM)920*920*93G 

886 930 IFCFQUD)931 >931 >932 

887 920 IFCFQUD)933*933*934 

888 933 RET8=SMLW 

889 60 T0 910 

891 934 RET 8=0 • 

892 60 T091O 

893 931 XI=FM*pI/WDTH 

894 RET8=CSINCXI*CFMU0+1 .0 ) +SML.W) -SINCXI*F MljD*SMLw) ) /XI 

896 60 T0 910 

897 932 JK=FQUD 

898 XA= ( - 1 «0)**JK 

899 XB=FM*p I+SMLW/WDTH 

901 XC=X3*FMUD 

902 XD=XB*CFMUD+1 .0) 

903 XE=FM*pI/WDTH 

904 XF=FQUD*PI/SMLW 

906 IF(XE~XF)960* 9 61*960 

907 961 RET8=SMLW*C0S<FQUD*FMUD*Pl)/2 .0 

908 60 T0 910 

909 960 XG=XE/(XF*XF-X£*XE> 

91 1 RET8=X6*CSIN<XC)-XA*SINCXD) ) 

912 910 CONTINUE 

913 firsum=firsum-epm*smusum*ret8/(ret*wdth> 

914 200 CONTINUE 

915 955 F0RMATC4F1O. 5*213) 

916 XI (QU2*MU2*MBAR)=FIRSUM 

918 176 CONTINUE 

919 175 CONTINUE 

920 IFCNN-NS)6*7*7 

922 6 00 188 «U2=1*MUMAX 

923 00 182 QU2=1*QUMAX 

924 RET9<QU2*MU2*MBAR)=X1 <QU2* MU2* MBAR ) 

925 182 CONTINUE 

926 188 CONTINUE 

927 NN=NN+ 1 

928 G0 TOl 

929 7 DO 183 MU2=l*MlJMAX 

930 D0 184 QU2=1*QUMAX 

931 A=A 

933 184 CONTINUE 

934 183 CONTINUE 

935 980 FORMAT <2HN0) 

938 956 F0RMATCF1O. 5*213) 

940 RETURN 
950 END 
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SSVEL2 


13:29 ROCKET 


50C PRESSURE D I ST R I BUT 1 0N-C0MP ART WENT SIDE 

55 INTEGER QU #QQ# QUMAX# QU2 # Q 1 MAX 

56 D I MENS I0N RET9C20#20# 1 ) # XI C20#20# 1 > 

70 10 PRINT 915 


72 915 F0RMATC 1 H0#2X#26HINPUT: 

73 INPUT #PHE£1 # ICON# NUMB* STEP 

75 PR I NT 916 

76 916 F0RMATC1H0#2X# 1 9HINPUT: 

77 INPUT#MU1# NS#Q1MAX 

78 PRINT 917 

79 917 F 0RMAT C 1 HO #2X#20HINpUT * 

80 INPUT# Y# STEP# SMLEL 1 

81 PRINT 918 

82 918 FORMAT C 1H0#2X# 5 1 HI Np UT : 
83&SMLW# SMLEL# XLNGTH) 

84 INPUT# MIJMAX# QUMAX# MMAX# MBAR 

87 PHEE=PHEE1/WDTH 

88 PI=3* 141592 

89 15 D0 170 MU2=1#MUMAX 

90 D0 171 QU2 = 1 # QUMAX 

92 FMBAR=FL0AT(M8AR) 

93 PI =3. 141 592 


PHEE# ICON# NUMB# STEP ) 

MUl # NS# Q 1 MAX ) 

Y# STEP# SMLEL1 ) 

MUMAX# QUMAX# MMAX# MBAR# WDtR# 
WDTH#SMLW# SMLEL# XLNGTH 


94 RET 9 C QU2# MU2# MBAR ) =0 • 

95 171 CONTINUE 

96 170 CONTINUE 

97 DO 172 MU2=1#M(JMAX 

98 0U2=1 

99 FMU2=FL0AT(MU2-1 ) 

101 FQU2=FL0AT CQU2- 1 ) 

102 IFCFMBARJ21 #21 #22 

103 22 IFCFQU2J31 #31 #32 

104 21 IF (FQU2>33#33#34 

105 33 RET9CGU2#MU2#MBAR>=SMLW 

107 GOTO 11 

108 34 RET9(QU2#MU2#MBAR)=0.0 

109 GO TO 11 

111 31 XI =FMBAR*p I / WDTH 

1 12 RET9CQU2#M(J2#MBAR) = CSIN<XI*<FMU2+1 »0)*SMLW) 

1 13&-SINCXI*FMU2*SMLW))/XI 

114 GO T0 11 

115 32 JK=FQU2 

116 XA= < - 1 • 0 > ** JK 

117 XF=FQU2*PI/SMLW 

118 XB=FM8AR*PI*SMLW/WDTH 

119 XC=XB*FMU2 

121 XD=XB*CFMU2+1 *0> 

122 XE=FMBAR*pi/WDTH 

123 IFCXE-XF)58# 59# 58 

124 59 RET9CQU2#MU2#MBAR)=SMLW*C0SCFMU2*FQU2*PI ) /2.0 
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125 G0 T0 11 

126 58 XG=X£/CXF*XF-XE*XE) 

128 RET9CQU2* MU2* M8AR ) =XG* < 3 I N< XC > -XA*S I NC XD > ) 

129 11 CONTINUE 

131 173 CONTINUE 

132 172 CONTINUE 

133 CALL ARICPHEE*WDTH*SMLEL*SMLW*XLNGTH,RET9* 

1 34&MUMAX*QUMAX * MMAX* MBAR* MUl * NS* Q 1 MAX* FMUD* F QUD* X 1 ) 

135 PI=3. 141592 

136 WDTH=20. 

140 MU=MU1 

150 FMU=FL0ATCMU-1 ) 

152 DUM=CFNU+1 »0)*SMLW 
160 19 QUSUM=0. 

170 D0 210 QU=1*QUMAX 
180 FQU=FL0ATCQU-1) 

190 IF C FQ U > 40 * 4 1 * 40 

200 41 EPQU=1»0 

210 60 T0 50 

220 40 EPQU =2 »0 

230 50 AB=FQU*PI/SMLW 

240 DISCR=PHEE*pHEE-A8*AB 

250 IF <DISCR>42* 43*44 

259 43 STOP 

260 42 AR1=SQRTC-DISCR)*SML.EL1 

261 AR4=SQRTC-DISCR>*SMLEL 

262 AR2=(EXPCAR1 )-EXP(-ARl ) >/2.0 

263 AR3=(EXPCAR4)+EXP(-AR4) )/2.0 
270 AT=-SQRTC-DISCR)*AR2 

280 AU=AR3 
290 60 T0 60 

300 44 at=sqrtcdiscr>*sin<sqrt<discr>*smleli ) 

310 AU=C0S(SQRT(DISCR)*SMLEL) 

319 FQp=0.0 

320 60 CALL ARGCFQP* FMU* FM8AR* SMLW* WDTH* RET > 

330 AV=C0SCFQU*PI*<Y“FMU*SMLW>/SMLW> 

340 qui=epqu*xi <QU>MU*MBAR)*AV/SMLW 

350 qusum=qusum+qui 

360 210 CONTINUE 

375 PRINT 900* Y * QUSUM 

380 IF<Y-DUM>70*71*71 

390 71 Y=DUM 

400 G0 T0 200 

410 70 Y=Y+.l 

41 1 G0 T0 1 9 

430 200 CONTINUE 

440 900 F 0RMAT C 1 H0*2X* 2HY = * F6»3*3X* 4HpBS=* F~ 1 0 • 4) 

450 GO T0 10 
460 END 

710 SUBROUTINE TANHY (X* Y* Z > 

715 RR=SQRT(A8S<X)) 
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720 RS*RR*Y 

725 Z*-RR*CEXPCRS>-EXP<-RS> >/CEXP<RS> + EXP<-RS> > 

730 RETURN 
735 END 

740 SUBROUTINE TAN<X*Y*Z> 

741 RT=SQRTCX> 

742 RU=RT*Y 

743 Z=RT*CSINCRU)/C0SCRU>> 

745 RETURN 

746 END 

747 SUBROUTINE ARG< TA* TB* TC * TD* TE* TF ) 

748 PI=3. 141592 

749 IFCTC)720*720*730 

751 730 IF<TA)731 *731 *7 32 

752 720 IFCTA1733* 733*734 

753 733 TF=TD 

754 GO TO 710 

755 734 TF=0. 

757 GO TO 710 

758 731 XI=TC*PI/TE 

7 59 TF = CSINCXI*<TB+1 .0 >*T‘D> -SINCXI*TB*TD) ) /XI 

760 GO TO 710 

761 732 JK=TA 

762 XA=C- 1 . 0 ) ** JK 

763 XB=TC*PI*TD/TE 

764 XC=XB*TB 

765 XD=XB*(TB+1 .0) 

766 XE=TC*PI/TE 

767 XF=TA*PI/TD 

768 IFCXE-XFJ7 60* 761*7 60 

770 761 TF=TD*G0S(TA*TB*PI)/2.O 

772 GO TO 710 

773 760 XG=XE/CXF*XF-XE*XE) 

774 TF=XG*(SIN(XC)-XA*SIN(XD) > 

776 710 CONTINUE 

777 RETURN 

778 END 

779 SUBROUTINE ARI <PHEE* WDTH* SMLEL* SMLW* XLNGTH* RET9* 
7804MUMAX *QUMAX * MMAX* MBAR* MU 1 * NS* Q1 MAX* FMUD* FQUD* X 1 > 

781 INTEGER QU*QQ*QUMAX*QU2*Q1MAX 

782 DIMENSION RET9C20*20* 1 ) *X1 (20*20* 1 > 

783 NN=1 

784 1 DO 175 MU2=1*MUMA£ 

785 DO 176 QU2=1*QUMAX 

786 Pl=3. 141592 

787 FIRSUM=0. 

7 88 DO 200 M= 1 * MMAX 

790 FM=FL0ATCM-1 ) 

791 AA=FM*pI/WDTH 

792 DISCR=PHEE*PHEE-AA*AA 

793 IF(DISCR)50* 51*52 
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794 50 ART=SQRTC-0ISCR3*XLNGTH 

796 IF CART-4 »0 > 850* 85 1 *851 

797 850 RR=SQKTCAriS<DISCR3> 

798 RS=RR*XLN6TH 

799 RET = -RR*CEXPCRS3-EXP<-RS3 )/CEXP<Ki>3 + EXPC-RS3 3 

801 60 T0 852 

802 851 RET = -SQRTCAriS<DlSCR> 3 

803 8 52 60 T03OO 

804 51 ST0P 

806 60 T0 10 

807 52 RT=5QRTC0ISCR) 

808 RU=RT*XLNGTH 

809 RET=RT*<SIN<RU 3/C0SCRU3 3 

811 300 SMUSUM=0 • 

812 00 220 MU S 1 * M.UMAX 

813 FMU=FL0AT<Mij-l 3 

814 QSUM=0. 

815 D0 210 Qll = 1 / Q 1 MAX 

816 FQU=FL0AT <QU 3 

817 FQU=FQU-1 *0 

818 I FCFQU360/ 60/ 62 

819 60 EPQU=1»0 

821 60 T0 301 

822 62 EP Q U =2 • 

823 301 A6=FQU*PI/SMLW 

824 DISCR2=PHEE*PHEE-A8*AB 

826 IF (DISCR2 365/66/67 

827 65 ARP = SORT ( -DI SCR2 3*SMLEL 

828 I F < ARP-4 •03860*861*861 

829 860 RR=SQRT(ABS<DISCR23 3 

830 RS=RR* SMLEL 

832 RET2 = -RR*(EXPCRS3-EXP<-RS3 3 /< EXP <KS 3 + EXP < -RS 3 3 

833 60 T0 347 

834 861 RET2=-SQRTCA8S(OISCR233 

836 60 T 0 347 

837 66 PRINT 980 

838 10 ST0P 

839 67 RT=SQKT C0ISCR2 3 

841 RU=RT +SMLEL 

842 R£T2=RT*(SIN(RU3/C08(RU3 3 

843 347 IF C Ft v J 3 720 / 720/ 730 

844 730 IFCFQU3731/731/732 

846 720 IFCFQU3733/ 733/734 

847 733 RET3=SMLW 

848 60 T0 710 

849 734 RET3=0 • 

851 60 T0 710 

852 731 XI=FM*PI/WDTH 

853 RET3=(SINCXI*(FMIJ+1 .0 3 +SMLW3 -SI NC XI *FM(J*SML W3 ) /XI 
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854 G0 T0 710 

856 732 JK=FQU 

857 XA=<-1 

858 XB=FM*PI*SMLW/WDTH 

859 XC=XB*FMU 

861 XD=XB*CFMU+1 • ) 

862 XE=FM*PI/WDTH 

863 XF = FQ(J*PI/SMLW 

864 IFCXE-XF >760* 7 6 1 * 7 60 

866 761 HET3=SMLw*C0S(FMU*FQU*PI )/2-0 

867 G0 T0 710 

868 760 XG=XE/<XF*XF-XE*XE) 

869 RET3=XG*CSIN<XC>-XA*SINCXD) > 

870 710 CONTINUE 

$71 FQUD=FL0AT<QU2-1 ) 

872 FMUO=FL0ATCMU2-n 

873 FMBAR=FL0ATCMBAR) 

874 QSUM=QSUM+EPQU*RET2*RET3*R£T9CQU>MU>MBAR) /SMLW 

876 210 CONTINUE 

877 SMUSUM=SMUSUM+QSUM 

878 220 CONTINUE 

879 IF C FM >70*70*71 

881 70 EPM= 1 • 

882 G0 T0 303 

883 71 EPM=2.0 

884 303 IF CFM) 920 *920 *930 

886 930 IF<FQUD>931 *931 *932 

887 920 IF CFQUD>933* 933* 934 

888 933 RET8=SMLW 

889 G0 T0 910 
891 934 RET8=0 • 

892. G0 T091O 

893 931 XI=FM*PI/WDTH 

894 RET8=<SINCXI*CFMU0+1 .0 >*SMLW> -SIN<XI*FMuD*SMLW> ) /XI 

896 G0 T0 910 

897 932 JK=FQUD 

898 XA=C-1 .0>**JK 

899 XB=FM*PI*SMLW/WDTH 

901 XC=XB*FMUD 

902 XD=XB*CFMUD+1 .Q> 

903 XE=FM*PI/WDTH 

904 XF=FQUD*Pl/SMLW 

906 IF (XE-XF >960*961*960 

907 961 RET8=SMLW*C0S<FQUD*FMuD*PI>/2 .0 

908 G0 T0 910 

909 960 XG=XE/CXF*XF-XE*XE> 

911 ret8=xg*csincxc>-xa*sin<xd> > 

912 910 CONTINUE 
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913 FIRSUM=FIRSUM-EPM*SMUSUM*RET8/(RET*WDTH> 

914 200 CONTINUE 

915 955 FORMAT C4F 10»5#2I3) 

916 XI CQU2>MU2^ MBAR)=FIRSUM 

918 176 CONTINUE 

919 175 CONTINUE 

920 IF CNN -NS )6*7 f 7 

922 6 DO 188 MU2=1>MUMAX 

923 DO 182 QU2= 1 * QUMAX 

924 RET9CQU2*MU2> MBAR)=X1 C Q U2 > Mij2 > MBAR ) 

925 182 CONTINUE 

926 188 CONTINUE 

927 NN=NN+ 1 

928 GO TO 1 

929 7 DO 183 MU2=1 /MUMAX 

930 DO 184 QU2=1>QUMAX 

931 A=A 

933 184 CONTINUE 

934 ,183 CONTINUE 

935 980 FORMAT C2HN0) 

938 956 FORMATCFIO.5^213) 

940 RETURN 
95t> END 
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SS8PD2 


11:28 ROCKET 


50 c pressure distribution-chamber side 

55 INTEGER QU>QQ> QUMAX> GU2 > Q 1 MAX 

56 DIMENSION RET9 (20>20> 1 ) > XI <20> 20> 1 > 

60 10 PRINT 915 

61 Y =0 • 

65 91 5 F0RMATC1H0>2X>26HINPUT: PH EE* ICON* NUMB* STEP > 

70 INPUT>PHEE1 > ICON* NUM8> STEP 

71 PRINT 917 

72 917 F0RMATC1HO>2X> 19HINPUT: MU1> NS> Q 1 MAX ) 

73 I NPUT> MU 1 > NS* Q 1 MAX 

75 IF ( I CON ) 20 > 20 > 30 

76 30 PRINT 916 

77 916 FORMATC 1H0>2X> 51HINpUT:MUMAX>QUMAX>MMAX> MBAR> WDTH> 
7 8&SMLW>SMLEL> XLNGTH ) 

79 input>mumax>qumax> mmax>mbar> wdth> sml w> smlel> xlngth 

81 20 PRINT 910 

82 N= 1 

83 PRINT 91 1 

84 PRINT 9 1 8* PHEE 1 > WDTH * SMLW> SMLEL> XLNGTH 

85 PRINT 912 

86 PRINT 9 1 9> MMAX> MUMAX> QUMAX> MBAR 

87 PHEE=PHEE1/WDTH 

89 15 DO 170 MU2= 1 > MUMAX 

90 DO 171 QU2 = 1 > QUMAX 

92 FMBAR=FLOAT (MBAR ) 

93 P I =3 • 1 4 1 592 

94 RET9 (QU2> MU2> MBAR ) =0 • 

95 171 CONTINUE 

96 170 CONTINUE 

97 DO 172 MU2= 1 > MUMAX 

98 QU2= 1 

99 FMU2=FL0AT(MU2-1 ) 

101 FQU2=FL0AT(QU2-1 ) 

102 IF CFMBAR >21*21*22 

103 22 IF(FQU2>31 >31 >32 

104 21 IFCFQU2)33>33>34 

105 33 RET9CQU2>MU2>MBAR)=SMLW 

107 GO TO 11 

108 34 RET9 <QU2> MU2 > MBAR ) =0 »0 

109 GO TO 11 

110 PRINT 911 

111 31 XI =FMBAR*P I /WDTH 

1 12 RET9<QU2>MU2>MBAR) = CSINCXI*CFM(J2+1 *0)*SMLW) 

1 13&-SIN(XI*FMU2*SMLW))/XI 

114 GO TO 11 

115 32 JK=FQU2 

116 XA=(-1 • 0 ) **JK 

117 XF=FQU2*PI/SMLW 

118 XB=fMBAR*Pl*SMLW/WDTH 
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119 XOXB+FMU2 

120 PRINT 912 

121 XD=XB*CFM(J2+1 .0) 

122 XE=FMBAR*PI/WDTH 

123 IFCXE-XF)58>59>58 

124 59 RET9CQU2>MU2>M6AR)=SMLW*C0S<FMU2*FQU2*Pl > /2»0 

125 G0 T0 11 

126 58 XG=XE/<XF*XF-XE*XE> 

128 RET9CQU2>MU2>MBAR)=XG*<SIN<XC)-XA*SIN(XD) > 

129 11 C0NTINUE 

131 173 CONTINUE 

132 172 CONTINUE 

133 CALL AR I (PHEE> WDTH> SMLEL> SMLW> XLNGTH> RET9 > 
134&MUMAX>QUMAX>MMAX>MBAR>MU1 >ns>qimax>fmud>fqud>xi > 
135 PI=3. 141592 

140 FIRSUM=0. 

145 D0 200 M=1>MMAX 
150 FM=FL0AT<M-1 ) 

155 AA=FM*PI/WDTH 

160 DlSCR=PHEE*pHEE-AA*AA 

165 IFCDISCR>50>51 >52 

170 50 ART=SQRT C -D I SCR ) *XLNGTH 

175 IFCART-4. )850>851 >851 

180 850 CALL TANHYCDI SCR> XLNGTH> RET > 

185 GO TO 852 

190 851 RET=-SQRTCABS<DISCR) ) 

195 852 GO TO 300 
200 51 PRINT 980 
205 GO TO 10 

210 52 CALL TAN<DISCR> XLNGTH>RET > 

215 300 QSUM=0. 

220 DO 210 QU= 1 > QUMAX 
225 FQU=FL0AT(QU-1 > 

230 IF(FQU)60>60>62 

235 60 EPQU=1 .0 

240 GO TO 301 

245 62 EpQU=2.0 

250 301 AB=FQU*PI/SMLW 

255 DISCR2=PHEE*PHEE-AB*AB 

260 IF CD ISCR2 )65>66>67 

265 65 ARP=SQRTC-DISCR2 )*SMLEL 

270 IF (ARP -4 •)860>861>861 

275 860 CALL TANHY ( DI SCR2 > SMLEL> RET2 ) 

280 GO TO 347 

285 861 RET2'’=-SQRT<ABSCDISCR2>> 

286 GO TO 347 
290 66 PRINT 980 
295 GO TO 10 

300 67 CALL TANCDISCR2> SMLEL> RET2 > 

305 347 SMUSUM=0»0 

310 302 DO 220 MU=1>MUMAX 
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315 FMU=FL0AT(MU-1 > 

325 CALL ARG(FQU*FMU*FM*SMLW* WDTH*RET3> 

330 FMBAR=FL0AT (MBAR ) 

334 FQP=0 »0 

343 SMUSUM=SMUSUM+EPQU*RET2*RET3*X1 CQU* MU* MBAR ) /SMLW 

345 220 CONTINUE 

355 QSUM=QSUM+SMUSUM 

360 210 CONTINUE 

365 IFCFM>70*70*71 

370 70 EPM= 1 • 0 

375 GO TO 303 

380 71 EPM=2.0 

385 303 FIRSUM=FIRSUM + EPM*QSUM*C0SCFM*PI*Y/WDTH) /CRET*WDTH) 

388 200 CONTINUE 

392 PRINT 994*PHEE1 *Y*FIRSUM 

393 994 F0RMATUH0>2X*5HPHEE=*F10.5>3X*2HY=*F5.2*3X*3HPC=*F12.5) 

394 IFCY-10. >90*91 *91 

395 90 Y=Y+ • 1 

396 GO TO 15 

397 91 GO TO 10 

660 910 F0RMAT<35HINITIAL VALUES FOR THIS CALCULATION//) 

665 911 F0RMATC1HO*3X*4HPHEE*4X*4HWDTH*4X*5HSMLW *3X*5HSML£L* 
666&2X* 6HXLNGTH) 

670 918 FORMAT CF8 • 4* 4F8 »2//) 

675 912 FORMAT C 1 HO* 4X* 4HMMAX* 3X* 5HMUMAX* 3X* 5HQUMAX* 4X* 4HM8AR ) 

680 919 F0RMAT<4C5X* I3>///) 

690 980 FORMAT C 9H0 1 SCR=0 *0 ) 

700 952 F0RMAT<I3*2C1PE12.4) ) 

705 END 

710 SUBROUTINE TANHYCX*Y*Z> 

715 RR=SQRT<ABS(X) ) 

720 RS=RR*Y 

725 Z=-RR*CEXP(RS)“EXP(-RS>)/<EXPCRS)+EXP(-RS) > 

730 RETURN 
735 END 

740 SUBROUTINE TAN(X*Y*Z> 

741 RT=SQRT<X) 

742 RU=RT*Y 

743 Z=RT*CSIN<RU>/C0SCRU>> 

745 RETURN 

746 END 

747 SUBROUTINE ARG( TA* TB* TC* TD* TE* TF ) 

748 PI=3- 141592 

749 IFCTO720* 720* 730 

751 730 IF(TA)731*731 *732 

752 720 IFCTA>733*733*734 

753 733 TF=TD 

754 GO TO 710 

755 734 TF=0. 

757 GO TO 710 

758 731 XI=TC*P I/TE 
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7 59 TF=CSINCXI*CTB+1 . 0 > *TD ) -S INC Xl+TB+TD) >/Xl 

760 G0 T0 710 

761 732 JK = TA 

762 XA=<-1 »0)**JK 

763 XB=TC*P I *TD/TE 

764 XC=XB*TB 

765 XD=XB*CTB+1 .0 > 

766 XE=TC*PI/TE 

767 XF=TA*pI/TD 

768 I F C XE-XF ) 76Q.» 7 61*760 

770 761 TF=TD*C0SCTA*TB*pi)/2.O 

772 60 T0 710 

773 760 XG=XE/CXF*XF-XE*XE) 

774 TF=XG*CSItM<XC>-XA*SlNCXD) ) 

776 710 CONTINUE 

777 RETURN 

778 END 

779 SUBROUTINE AR I CPHEE* WDTH> SMLEL* SMLW* XLNGTH* RET9* 

780&mumax*qumax*mmax*mbar.,mui >ns*qimax*fmud*fqud*xi ) 

781 INTEGER QU*QQ*QUMAX*QU2*Q1MAX 

782 DIMENSION RET9C20*20* 1 ) * XI C20*20> 1 ) 

783 NN= 1 

784 1 DO 175 MU2= 1 * MUMAX 

785 DO 176 QU2=1,QUMAX 

786 Pl=3* 141592 

787 FIRSUM=0. 

788 DO 200 M= 1 » MMAX 

790 FM=FL0AT CM- 1 ) 

791 AA=FM*P I / WDTH 

792 D I SCR=PHEE*PHEE-AA*AA 

793 IFCDISCR) 50 >51^52 

794 50 ART=SQRT<-DISCR>*XLNGTH 

796 IF < ART-4 »0 ) 850> 85 1 >851 

797 850 RR = SQRT<ABS<DISCR) ) 

798 RS=RR*XLNGTH 

799 RET=-RR*CEXP<RS)-EXP<-RS) )/CEXP(RS)+EXPC-RS) > 

801 GO TO 852 

802 851 RET = ~SQRT CABS CDI SCR ) ) 

803 8 52 GO T0300 

804 51 STOP 

806 GO T 0 10 

807 52 RT=SQRTCDISCR> 

808 RU=RT*XLNGTH 

809 RET=RT*CSIN<RU)/C0SCRU> > 

811 300 SMUSUM=0. 

812 D0 220 MU=1>MUMAX 

813 FMU=FL0ATCMU-1 3 

814 QSUM=0. 

815 DO 210 QU = 1 •» Q 1 MAX 

816 FQU=FL0ATCQU) 

817 FQU=FQU-1 .0 
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818 IFCFQU>60>60^62 

819 60 EPQU=1 .0 

821 G0 T 0 301 

822 62 EPQU=2. 

823 301 AB=FQU*PI/SMLW 

824 DISCR2=PHEE*PHEE-AB*A8 

826 IF (DISCR2 >65.»66*67 

827 65 ARP = SORT (-DISCR2 >*SMLEL 

828 IFCARP-4. 0)860>861 ^861 

829 860 RR = SQRT CABS (DISCR2 ) ) 

830 RS=RR*SMLEL 

832 RET2=-RR*CEXPCRS)-£XPC-RS> > /C EXP CRS >+ EXP C -RS > > 

833 G0 T 0 347 

834 861 RET2=-SQRTCABSCDISCR2) ) 

836 G0 T0 347 

837 66 PRINT 980 

838 10 ST0P 

839 67 RT=SQRTCDISCR2) 

841 RU=RT*SMLEL 

842 RET2=RT*(SIN(RU>/C0SCRU> > 

843 347 IFCFM)720> 720^730 

844 730 IF(FQU>731-*731.» 732 

846 720 IF CFQU)733> 733^734 

847 733 RET3=SMLW 

848 G0 T 0 710 

849 734 RET3=G. 

851 G0 T0 710 

852 731 X I =FM*p I /WDTH 

853 RET3=CSIN(XI*CFMU+1 *0 >*SMLW>-SINCXI*FMij*SMLW> > /XI 

854 G0 T0 710 

856 732 JK=FQU 

857 XA=C- 1 .05**JK 

858 XB=FM*Pl*S«LW/WDTh 

859 XC=X8*FMU 

861 XD=X8*CFMU+1 . > 

862 XE=FM*pI/WDTH 

863 XF=FQU*PI/SMLW 

864 IF ( XE-XF ) 7 60.* 76 1 j 7 60 

866 761 RET3=SMLW*C0S<FMU*FQU*PI>/2.O 

867 60 T0 710 

868 760 XG=XE/CXF*XF-XE*XE> 

869 RET3=Xe*CSIN(XC>-XA*SIN(XD> > 

870 710 CONTINUE 

871 FQUD=FL0AT(QU2-1 ) 

872 FMUD=FL0ATCMU2-1 ) 

873 FMBAR=FL0AT CMBAR ) 

874 QSUM=QSUM+EPQU*RET2*RET3*RET9(QU» MU^MBAR) /SMLW 

876 210 CONTINUE 

877 SMUSUM=SMUSUM+QSUM 

878 220 CONTINUE 

879 IF (FH>70* 70*7 1 
881 70 EPM=1 . 
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882 G0 T 0 303 

883 71 EPM=2 • 0 

884 303 IF(FM) 920*920*9 30 

886 930 IFCFQUD)931 >931 >932 

887 920 IF (FQUD) 933*933*934 

888 933 RET8=SMLW 

889 G0 T0 910 

891 934 RET8=0. 

892 G0 T091O 

893 931 Xl=FM*pI/WOTH 

894 RET8=CSlNOa*CFMUD+l .0 ) *SML W) -S IN CX I *FMUD*SML W> ) /X I 

896 G0 T0 910 

897 932 J«=FQUD 

898 XA=<- 1 • 0 ) ** JK 

899 XB=FM*PI*SMLW/WDTH 

901 n XC=XB*FMUD 

902 XD=XB*CFMUD+1 .0) 

903 XE=FM*Pl/WDTh 

904 XF=FQUD*PI/SMLW 

906 IFCXE~XF)960* 961 *96Q 

907 961 RET8=SMLW*C0SCFQUD*FMUD*PI)/2.O 

908 G 0 T0 910 

909 960 XG=XE/CXF*XF-XE*XE) 

911 RET8=XG*<SIN(XC)-XA*SINCXD) ) 

912 910 CONTINUE 

913 FIRSUM=FIRSUM-£PM*SMUSUM*RET8/CRET*WDTH) 

914 200 CONTINUE 

915 955 F0RMATC4F1O. 5*213) 

916 XI CQU2*MU2*MBAR)=FIRSUM 

918 176 CONTINUE 

919 175 CONTINUE 

920 I F ( NN-NS >6* 7*7 

922 6 DO 188 MU2=1*«UMAX 

923 DO 182 QU2=1*QUMAX 

924 RET9(QU2*MU2*M8AR)=X1 < QU2 * MU2 * MBAK ) 

925 182 CONTINUE 

926 188 CONTINUE 

927 NN=NN+ 1 

928 GO T01 

929 7 DO 183 Mu2=l*tMUMAX 

930 DO 184 QU2=1*QUMAX 

931 A=A 

933 184 CONTINUE 

934 183 CONTINUE 

935 980 F 0RMAT C2HN0 ) 

938 956 F0RMATCF10. 5*213) 

940 RETURN 
950 END 
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11:09 ROCKET 


50 C CHARACTERISTIC EQ’tM FOR THE BAFFLE PROBLEM CACTI VE BOUNDARY) 

51 COMPLEX R£T9#DISCR#PHEE#ART#RET#SMUSUM#QSUM#DISCR2#ARP#RET2# 
52&DSUM#FISUQU#SECSUM#DISCR3#ARS#RET5#X1#T0SUM#PHEE2#PHEE3#FIRSUM# 
53&PHEE1 #ADA#REA#REB# YMA# YQ#R£C#R£D#RE£#REF#FIND #DFEFUN#YQ1 # YQ2 
55 I NTEGER QU#QQ# QUMAX# QU2# Q1 MAX 

60 DIMENSION RET9C20#20# 1 > # XI <20#20# 1 > #T0SUMC 3) 

65 10 PRINT 915 

67 915 FORMAT C 1H0#2X#26HINPUT : PHEE# I CON, NUMB, STEP ) 

69 INPUT#PHEE1 # I CON# NUMB# STEP 
71 IFCIC0N>20#20#30 
73 30 PRINT 916 

75 916 F0RMATC1H0#2X# 1 9HINPUT: MU 1 # NS# Q 1 MAX ) 

77 I NpUT # MU 1 # NS#Q 1 MAX 
79 PRINT 917 

81 917 FORMAT C 1 H0#2X# 1 3HINPUT : YQ#YMA) 

83 INPUT# YQ# YMA 
85 PRINT 918 

87 918 F ORMAT C 1 HO #2 X# 51 HI NPUT : MyMAX# QUMAX# MMAX# MBAR# WDTH# SMLW# 
88&SMLEL# XLNGTH) 

90 INPUT#MUMAX# Q UMAX # MMAX# MBAR# WDTH# SMLW# SMLEL# XLN GTH 
110 PHEE=PHEE1 /WDTH 

120 20 N= 1 

121 PRINT 913 

146 15 DO 169 J J= 1 # 3 

147 F J J=FL0AT C J J - 1 ) 

148 FIND=CMpLXC .000001 # .000001 > 

149 YQ=YQ+FJJ*FIND 

155 DO 170 MU2=1*MUMAX 
160 D0 171 QU2= 1 # QUMAX 
165 FMBAR=FL0ATCMBAR) 

170 PI=3. 141592 

175 RET9CQU2#MU2#MBAR)=CMpLXC0.#0. ) 

180 171 CONTINUE 
185 170 CONTINUE 
190 DO 172 MU2=1#MUMAX 
195 QU2=1 

200 FMU2=FL0ATCMU2-1 ) 

205 FQU2= FLOAT (QU2- 1 ) 

210 IFCFMBAR)21 #21 #22 

215 22 IF C FQU2 )31#31#32 

220 21 IF CFQU2>33# 33# 34 

225 33 RET9CQU2#MU2#MBAR)=SMLW 

230 GO TO 11 

235 34 RET9CQU2#MU2#MBAR)=0*0 

240 GO TO 11 

250 31 XI=FMBAR*PI/WDTH 

255 RET9CQ|J2#MU2#MBAR) = CSINCXI*CFMU2+ 1 .0>*SMLW) 
260&-SlNCXl*FMU2*SMLW))/XI 
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265 G0 T0 11 

270 32 JK=FQU2 

275 XA=<-1 .0>**JK 

280 XF=FQU2*PI/SMLW 

285 XB=FMBAR*pl*SMLW/W0TH 

290 XC=XB*FM(J2 

300 XD=X8*CFMU2+1 .0) 

305 XE=FMBAR*pI/WDTH 
310 IFCXE~XF)58*59i 58 

315 59 RET9CQU2*MU2>M8AR)=SMLW*C0SCFMU2*FQU2*PI ) /2.0 
320 G0 T 0 11 

325 5B XG=XE/CXF*XF-XE*XE> 

330 RET9CQU2#MU2»MBAR)=XG+CSIN(XC)-XA*SINCXD) ) 

335 11 CONTINUE 
3.40 173 CONTINUE 
345 172 CONTINUE 

350 CALL ARICPHEE> WDTH>SMLEL>SMLW*XLN6TH#RET9>Y«A>YQv* 
355&MUMAXjQUMAX*MMAX>MBAR>MU1 »NS>Q1MAX> FMUD*FQUD.» XI > 

360 PI=3. 141592 

365 FIRSUM=CMPLXC0.>0.) 

370 D0 200 M= 1 > MMAX 
375 FM=FL0AT CM- 1 ) 

380 AA=FM*PI/WDTH 

385 D I SCR=PHEE*PHEE-AA*AA 

390 IF (CABS CD I SCR ) ) 50 » 51 » 52 

395 50 ART=CSQRTC-DISCR)*XLNGTH 

420 852 G0 T0 300 

425 51 PRINT 980 

430 G0 T0 10 

435 52 CALL TAN CD I SCR* XLNGTH* RET ) 

436 ADA=CMPLX(0 •» I •) 

437 REA=RET-ADA*PHEE*YMA 

438 REB=1 .0+ADA*pHEE*YMA*RET/(CSORTCDISCR)*CSQRTCDlSCR) > 

439 R£T=REA/REB 

440 300 SMUSUM=CMPLX(0.»0.) 

445 D0 220 MU=1 » MUMAX 

450 FMU=FL0AT<MU-1 > 

455 QSUM=CMPLX (0 • > 0 • ) 

465 00 210 QU= 1 » QUMAX 
470 FQU=FL0ATCQU-1 ? 

475 I F ( FQU 1 60* 60# 62 

480 60 EPQU=1 «0 

485 G0 T0 301 

490 62 EPQU =2 »0 

495 301 AB=FQU*PI/SMLW 

500 D1SCR2=PHEE*PK£E-A6*AS 

505 IFCCABSCDISCR2) )65*66>67 

510 65 ARP = CSQRT (-DISCR2 ) 

535 G0 T0 347 
540 66 PRINT 980 
545 G0 T0 10 


C-27 



550 67 CALL TANCDISCR2* Si v iLEL* RET2 ) 

551 REC=RET2-ADA*pHEE*YG 

552 RED= 1 .0+ADA*pHEE*YQ*RET2/CCSQRTCDISCR2>*CSQRTCDlSCR2> ) 

553 RET2=REC/RED 

555 347 F MBAR = FLOAT (MBAR ) 

560 CALL ARGCFGU>FMU^FM>SMLW^ WDTH>RET3> 

565 FQR=0 • 0 
570 FMU0=FMU 
575 FQUD=FQU 

580 QSUM=QSUM+EPQU*RET2*RET3*X1 C QU> MU* MBAR ) /SMLW 
585 

590 FMU1=FL0AT<MU1-1> 

595 210 CONTINUE 

600 IFCFMU-FMU1 >4*5>4 

605 5 DSU«=OSUW 

610 4 SMUSUM=SMUSUM+QSUM 

615 220 CONTINUE 

620 IFCFM)70>70,71 

625 70 EPM= 1 • 0 

630 GO TO 303 

635 71 EPM=2 • 0 

640 303 FIRSUM=FIRSUM+SMUSUM*SMUSUM*EPM/CRET*WDTH> 

645 200 CONTINUE 

650 SECSUM=CMPLX(0.>0. ) 

655 DO 230 MlJ = l>MU«AX 
665 FMtJ=FL0ATCMu-l ) 

670 F I SUGU = CMpLX C0*»0* ) 

675 DO 240 QU=1*QUMAX 
680 FQ,U=FL0ATCQU-1 ) 

685 IFCFQUi 80 >80^81 

690 80 EP Q U = 1 *0 

695 60 TO 500 

700 81 EPQU=2.0 

705 500 AC=FQU*PI/SMLW 

7 10 D I SCR3=PH£E*PHEE-AC*AC 

715 IF (CABS CD I SCR 3) )90>9 1 >92 

720 90 ARS=CSQRTC-DISCR3)*SMLEL 

750 91 PRINT 980 

755 GO TO 10 

760 92 CALL TAN(DISCR3* SMLEL* RETS) 

761 R£E=RET5-ADA*PHEE*YQ 

762 REF = 1 .0+ADA*pHEE*YQ*RET5/CCSQRTCDlSCR3)*CSQRTCDISCR3> > 

763 RET 5= REE/REF 
765 410 F GP=0 *0 
770 FMUD=FMU 
775 FQUD=FQU 

780 F ISUQU=FISUQU+RET5*X1 ( QU> MU-» MBAR ) *X 1 ( Q{J> MU* MBAR ) * 

785&EPQU/SMLW 

790 240 CONTINUE 

795 SECSUM=SECSUM+FISUGU 

796 230 CONTINUE 
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800 953 F0RMATCF12.6) 

805 T0SUM(JJ)=FIRSUM + SECSUM 

806 169 CONTINUE 

807 YQ2=YQ-FIND 

808 DFEFUN=(T0SUM(3)-T0S|JM(.1 ) >/C2.0*FIND> 

809 YQ=YO2-T0SUM(2)/DFEFUN 

810 YQ 1 =YQ 

815 PRINT 952*N*YQ1 *T0SUMC2 > 

816 DFEDUM=CA8S<T0SUMC2>> 

817 IFCDFEDUM- .000001 > 97 6* 976* 969 
820 969 I F ( N- NUMB ) 97 5* 975*976 
825 975 N=N+ 1 

830 YQ=YQ-FIND 
835 G0 T0 15 
840 976 G0 T0 10 

895 913 F 0RMAT <1H0*1X*1HI* 4X* 6HR-PHEE* 5X* 6H I -PHEE 
8964* 6X* 8HR-FEFUNC* 6X* 8HI -FEF UNC > 

910 980 F0RMATC9HDISCR=O*O) 

945 952 FORMAT <13* 4F 12*8) 

950 END 

985 SUBROUTINE TAN(X*Y*Z) 

986 COMPLEX X*RT*RU*Z 
990 RT=CSQRTCX) 

995 RU=RT*Y 

1000 Z=RT*CCSINCRU)/CC0SCRU) ) 

1005 RETURN 
1010 END 

1015 SUBROUTINE ARG CTA* TB* TC* TD* TE* TF ) 

1020 PI=3. 141592 

1025 IFCTC1720* 720* 730 

1030 730 IF(TA)731*731*732 

1035 720 IFCTA>733*733*734 

1040 733 TF=TD 

1045 GO T0 710 

1050 734 TF=0. 

1055 GO TO 710 
1060 731 XI=TC*PI/TE 

10.65 TF=(SINCXI*CT8+1 .0 ) *TD > -SI N<XI *TB*TD > ) /XI 

1070 GO TO 710 

1075 732 JK=TA 

1080 XA= C - 1 • 0 ) ** JK 

1085 XB=TC*PI*TD/TE 

1090 XC=XB*TB 

1095 XD=XB*<TS+1 .0) 

1100 XE=TC*PI/TE 

1105 XF=TA*PI/TD 

1110 I F (XE-XF >760*761*7 60 

1115 761 TF=TD*C0S<TA*TS*PI )/2 *0 

1120 GO TO 710 

1125 760 XG=XE/CXF*XF-XE*XE> 

1130 TF=XG*(SIN(XC)-XA*SIN<XD) ) 
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1135 710 CONTINUE 
1 140 RETURN 
1145 END 

1150 subroutine arkphee, wdth>smlel.,smlw,xlngth,ret9* yma,yq, 
1 1 55&MUMAX*QUMAX*MMAX* MBAR* MUl * NS* Q 1 MAX* FMUD* FQUD* XI > 

1156 COMPLEX FIRSUM,DISCR*ART*RR*RS*RET*RT*RU*SMUSUM.» 

1 1 57&QSUM* D ISCR2.» ARP* RET2 * RET9* X 1 .» ADA*REA* REB* YMA* YQ* R EC* RED 

1160 INTEGER QU*QQ*QUMAX*GU2*Q1 MAX 

1165 DIMENSION RET9C20*20* 1 ) * XI C20*20* 1 ) , 

1170 NN= 1 

1175 1 DO 175 MU2= 1 * MUMAX 
1180 DO 176 QU2=1*QUMAX 
1 185 pl=3. 141592 
1190 FIRSUM=CMPLXC0.*0. > 

1195 DO 200 M= 1 * MMAX 
1200 FM=FLOAT CM- 1 ) 

1205 AA=FM*PI/WDTH 
1210 DISCR=PHEE*PHEE-AA*AA 
1215 IF (GABS CD I SCR >) 50* 51 * 52 
1220 50 ART=CSQRTC-DISCR)*XLNGTH 
1250 851 RET = -CSQRT ( -D I SCR 3 
1255 852G0 T0300 
1260 51 STOP 
1265 GO TO 10 

1270 52 CALL TANCD ISCR* XLNGTH* RET ) 

1281 ADA=CMPLX<0.0* 1 .0> 

1282 REA=RET-ADA*pHEE*YMA 

1283 REB=1 .0+ADA*pHEE*YMA*RET/CCSQRT(DlSCR)*CSQRTCDISCR) > 

1284 RET=REA/REB 

1285 300 SMUSUM=CMPLXC0.#0. > 

1290 DO 220 M|J=1> MUMAX 

1295 FMU=FL0ATCMU-1 > 

1300 QSUM=CMPLX<0.*0. > 

1305 DO 210 QU=1*Q1MAX 
1310 fqu=floatcqu) 

1315 FQU=FQU-1 *0 
1320 IFCFQU160* 60*62 
1325 60 EPQU=1 *0 
1330 GO TO 301 
1335 62 EPQU=2. 

1340 301 A8=FQU*P I./SMLW 
1345 D1SCR2=PH£E*PHEE-AB*AB 
1350 IF (CABS (DISCR2 >>65*66*67 
1355 65 ARP=CSQRTC-DISCR2)*SMLEL 
1395 66 PRINT 980 
1400 10 STOP 

1405 67 CALL TANCDI SCR2* SMLEL*RET2 > 

1415 REC=RET2-ADA*PHEE*YQ 
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1416 RED*1 .0+ADA*pHEE*YQ*RET2/CCSQRT<DlSCrt2)*CSQRT<DlSCR2> > 

1417 RET2=REC/RED 

1420 347 CALL ARG< FQU* FMU-» FM* SMLW* WDTH* RET3 > 

1535 FQUD=FL0AT(QU2-1 ) 

1 540 FMUD=FL0AT (MU2- 1 > 

1545 FM8AR=FL0ATCM8AR) 

1 550 QSUM=QSUM+EPQU*RET2*RET3*RET9(QU>MU^MBAR) /SMLW 

1555 210 CONTINUE 

1560 SMUSUM=SMUSUM+QSUM 

1565 220 CONTINUE 

1570 IF CFM> 70* 70* 7 1 

1575 70 EPM=1. 

1580 GO TO 303 
1585 71 EPM=2.0 

1590 303 CALL ARG<FQUD*FMUD*FM* SMLW* WDTH*RET8) 

1600 

1705 FIRSUM=FIRSUM-EPM*SMUSUM*RET8/CRET*WDTH> 

1710 200 CONTINUE 

1720 XI <QU2>MU2#M8AR)=FIRSUM 

1725 176 CONTINUE 

1730 175 CONTINUE 

1735 IF C NN-NS "> 6 * 1*1 

1740 6 DO 188 MU2=1*MUMAX 

1745 DO 182 QU2= 1 * QUMAX 

1750 RET9 CQU2* MU2* MBAR ) =X 1 ( QU2* p)U2* MBAR ) 

1755 182 CONTINUE 
1760 188 CONTINUE 
1765 NN=NN+1 
1770 GO T01 
i 77 5 7 CONTINUE 
1800 980 FORMAT C2HN0) 

1810 RETURN 
1815 END 
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